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Abstract 

The  formation  of  cracks  in  elastic  media  such  as  sea  ice  generates  elastic 
waves  in  a  radiation  pattern  being  dependent  on  the  actual  fracture  process 
and  the  stratification  of  the  medium.  In  the  case  of  horizontal  stratification, 
this  phenomenon  can  be  idealized  and  mathematically  modeled  describing  the 
directionality  of  the  acoustic  emission  produced  by  compact  cracks  in  such  an 
environment. 

The  thrust  of  the  present  research  has  been  to  develop  an  analytical  and 
numerical  model  of  the  elastic  wave  field  in  range  independent  elastic  environ¬ 
ments  for  various  seismic  source  mechanisms.  The  source  types  being  considered 
are  explosive  sources,  point  forces,  shear  cracks,  and  tensile  cracks.  First,  the 
compact  source  representations  with  fault  surface  in  an  arbitrary  direction  will 
be  derived,  and  incorporated  in  a  numerical  model  for  propagation  in  strati¬ 
fied  elastic  media  to  yield  the  basic  Green’s  function  solution.  This  solution 
is  then  applied  to  derive  the  seismo-acoustic  field  produced  by  more  complete 
cracking  mechanisms,  like  non-compact  and  moving  cracks.  Finally,  the  effect  of 
anisotropy  on  the  acoustic  emission,  e.g.  in  sea  ice  and  periodically  fine  layered 
sediments,  will  be  considered. 

The  developed  model  is  applied  to  the  ice  cracking  noise  radiation  in  the 
central  Arctic  environment,  and  the  characteristics  of  the  field  produced  by  dif¬ 
ferent  source  and  environmental  parameters  are  discussed.  The  developed  model 
can  be  applied  to  the  source  inversion  problem,  i.e.  the  characterization  of  the 
cracking  mechanism  from  its  acoustic  emission  with  the  purpose  of  obtaining  a 
better  understanding  of  the  general  ambient  noise  in  the  central  Arctic.  Another 
expected  application  is  the  development  of  remote  sensing  techniques  suitable 
for  the  study  of  mechanical  properties  of  ice. 

Thesis  supervisor  :  Dr.  Henrik  Schmidt 
Associate  Professor  of  Ocean  Engineering 
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Glossary 


*  :  convolution. 

{ )  :  thickness-weighted  average. 

A  :  fault  area. 

A„  :  coefficient  matrix  for  n-th  layer. 

B™  :  unknown  variables  matrix  for  n-th  layer  and  m-th  Fourier  order. 

Cij  :  elastic  constants  for  a  transversely  isotropic  medium  (refer  to  Eq  5.2). 
Cp  :  phase  speed  of  the  P  wave  in  the  transversely  isotropic  medium. 
cSv  ’  phase  speed  of  the  SV  wave  in  the  transversely  isotropic  medium. 
csh  ’  phase  speed  of  the  SH  wave  in  the  transversely  isotropic  medium. 
cc  :  compressional  wave  velocity. 
ct  :  shear  wave  velocity. 

ei)  e2>  63  :  unit  vectors  in  (xx,  x2,  x3)  rectangular  coordinate  system. 

ez,  ev,  ez  :  unit  vectors  in  (x,  y,  z)  rectangular  coordinate  system. 

F m(r,z)  :  field  parameter  matrix  for  particular  solution  used  for  boundary 
condition  when  z  =  0 ,zn,  and  for  solution  at  receiver  depth  when  z  =  zr. 

Fm(s,z)  :  Hankel  transform  of  Fm(r, z). 

f  :  vector  body  force  used  in  the  representation  theorem. 

fn  :  force  in  n-direction. 

F  :  vector  body  force  as  a  source  term  in  the  equation  of  motion,  F  =  Fxex  + 
Fy&y  +  Fzez. 

F,  G,  H  :  displacement  potential  in  Cartesian  coordinates. 
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Gnp  :  Green’s  function  for  displacement  in  n-direction  due  to  a  body  force  in 
p-direction. 

GnPiq  :  derivative  of  Gnp  in  ^-direction. 

Hm(r,  z)  :  Inverse  Hankel  transform  of  Hm(s,  z),  and  m-th  order  field  param¬ 
eter  matrix. 

Hm(s,  z)  :  Hankel  integrand  modulus  matrix  for  a  set  of  field  parameters  used 
for  boundary  condition,  z  =  0,  zn  for  boundary  condition,  and  z  =  zr  for 
solution  at  receiver  depth. 

h  :  compressional  wavenumber. 

k  :  shear  wavenumber. 

L,  M ,  N  :  force  potentials  in  Cartesian  coordinates. 

Mo(t)  :  seismic  moment  source  function, 
n  :  normal  vector  (  usually  to  the  crack  surface  ). 
n,-  :  component  of  normal  vector  in  i-direction. 

Rm(r,  z)  :  Inverse  Hankel  transform  of  Rm(s,  z),  i.e.  m-th  order  field  param¬ 
eter  matrix  due  to  source. 

R m(s,z)  :  matrix  for  a  set  of  field  parameters  used  for  boundary  condition 
due  to  source,  z  =  0,  zn  for  boundary  condition,  and  z  —  zr  for  particular 
solution  at  receiver  depth. 

r  :  position  vector. 

r  :  range  in  rectangular  coordinates. 

s  :  horizontal  wave  number. 

[T]  :  traction  discontinuity. 

[u]  :  vector  displacement  discontinuity  on  the  crack  surface, 

u  :  vector  displacement. 

u  :  r-direction  displacement  at  horizontal  interface  in  cylindrical  coordinates. 
un  :  displacement  in  n  direction. 
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u(t)  :  average  slip  in  fault  source. 

v  :  0-direction  displacement  at  horizontal  interface  in  cylindrical  coordinates. 
w  :  ^-direction  displacement  at  horizontal  interface  in  cylindrical  coordinates. 
X0(£)  :  force  source  function. 

X,  Y,  Z  :  three  components  of  force  F  as  a  source  term  in  the  equation  of 
motion. 

zn  :  thickness  of  n-th  layer. 
zr  :  receiver  depth. 
za  :  source  depth. 

a  :  vertical  wavenumber  parameters  for  compressional  waves. 

{3  :  vertical  wavenumber  parameters  for  shear  waves. 

6  :  dip  angle,  refer  to  Fig  3.2. 
e  :  strain. 

A  :  rake  angle,  refer  to  Fig  3.2. 

A,  p.  :  Lam4  constants. 
p  :  density. 
a  :  stress. 

ozz  :  normal  stress  in  2-direction  at  horizontal  interface. 
arz  :  shear  stress  in  r-direction  at  horizontal  interface. 
oez  :  shear  stress  in  0-direction  at  horizontal  interface. 

4>s  :  strike  angle,  refer  to  Fig  3.2. 

4> ,  fa,  ^>2)  03  :  potentials  in  ( Xi,X2,Xs )  Cartesian  coordinate  system  used  in 
Chapter  2  and  3,  of  which  relation  to  displacement  are  defined  in  Eq  2.28. 

4>,  ipy,  4>z  '•  potentials  in  (x,y,2)  Cartesian  coordinate  system  used  in 

Chapter  3. 
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<j>,  Vv,  ipz  :  potentials  in  [r,6,z)  cylindrical  coordinate  system,  of  which 
relation  to  displacement  are  defined  in  Eq  3.11. 

<f>,  A,  V’  :  scalar  potentials  in  (r,  0,  z)  cylindrical  coordinate  system,  of  which 
relation  to  displacement  are  defined  in  Eq  3.1. 


£  :  parameter  defined  as 


f  =  sign(z  — 


> 


z  —  zt  >  0 
Z  -  Zg  <  0  . 
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Chapter  1 
Introduction 


1.1  Background 

The  formation  of  cracks  in  elastic  media  such  as  sea  ice  generates  elastic  waves 
in  a  radiation  pattern  being  dependent  on  the  actual  fracture  process  and  the 
stratification  of  the  medium.  In  the  case  of  horizontal  stratification,  this  phe¬ 
nomenon  can  be  idealized  and  mathematically  modeled  as  a  compact  but  di¬ 
rectional  source  radiating  in  a  stratified  elastic  medium,  yielding  the  Green’s 
function  solution.  Further,  the  basic  Green’s  function  solution,  by  integration 
in  time  and  space,  can  lead  to  synthesis  of  the  wave  field  produced  by  more 
complete  cracking  mechanisms,  like  non-compact  and  moving  cracks.  A  math¬ 
ematical  description  of  this  phenomenon  is  essential  to  the  inverse  problem  of 
determining  the  source  mechanism  from  experimental  data,  provided  that  the 
medium  is  known  on  the  basis  of  geological  and  geophysical  data,  and  the  for¬ 
ward  problem  of  predicting,  with  accumulated  data,  the  field  parameters  of 
interest  at  specific  positions.  Further,  such  a  mathematical  model  is  instrumen¬ 
tal  to  the  medium  inverse  problem  of  seismo-acoustic  exploration,  where,  for 
various  types  of  excitation,  the  observed  data  are  used  to  probe  the  geological 
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composition.  On  this  background,  the  research  described  in  this  thesis  is  cen¬ 
tered  around  the  development  of  an  efficient  mathematical/numerical  model  of 
the  acoustic  emission  produced  by  compact  and  propagating  cracks  in  realistic 
stratified  environments. 

To  illustrate  how  the  present  research  can  be  instrumental  in  understanding 
an  important  physical  problem,  an  application  to  the  central  Arctic  ambient 
will  be  discussed.  The  main  source  of  the  central  Arctic  ambient  noise  has  been 
known  to  originate  from  the  elastic  motion  of  the  ice  plate  caused  by  environ¬ 
mental  stresses.  When  the  stress  in  the  ice  plate  is  locally  greater  than  the 
strength  of  the  ice  causing  fracture,  the  stored  energy  is  released  in  the  form  of 
elastic  motion  of  the  ice.  This  energy  radiates  into  the  water,  forming  ambient 
noise  when  the  events  are  aggregated.  In  order  to  better  understand  the  Arctic 
ambient  noise  generating  mechanism,  two  levels  of  approaches  are  being  pursued 
[9].  One  approach  is  to  correlate  the  spectral  and  temporal  characteristics  of  am¬ 
bient  noise  to  gross  environmental  parameters  such  as  thermal  changes,  current 
and  wind  stresses  [33].  The  other  level  of  approach  is  to  look  into  the  individ¬ 
ual  ice  cracking  events,  and  treat  them  as  a  noise  generating  source  element, 
which  forms  the  average  ambient  noise  when  aggregated.  The  latter  approach  is 
based  on  the  assumption  that  the  central  Arctic  ambient  noise  is  dominated  by 
the  radiation  due  to  mechanical  processes  in  the  ice  cover,  so  that  the  ambient 
noise  is  aggregate  of  individual  events  [9]  [34].  Therefore,  it  is  obvious  that  the 
models  that  represent  probable  source  mechanisms  of  a  single  event  need  to  be 
developed  to  understand  the  central  Arctic  ambient  noise  characteristics  as  an 
aggregate  of  events. 

In  order  to  be  more  specific  about  the  issue  investigated  in  this  study  in  rela¬ 
tion  to  the  central  Arctic  ambient  noise,  the  processes  involved  in  the  generation 
of  the  Arctic  ambient  noise  can  be  categorized  as  (Fig  1.2) 
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(1)  Development  of  environmental  stresses, 

(2)  Ice  plate  motion  induced  by  fracture  due  to  the  environmental  stresses, 

(3)  Radiation  from  the  ice  plate  into  water. 

The  first  phenomenon  has  been  extensively  studied  and  associated  with  the  ob¬ 
served  ambient  noise  [31]  [33].  The  radiation  mechanism  from  the  ice  plate  to 
water  is  also  studied  by  some  authors  [10]  [42]  [24]  [23]  [22].  However,  previ¬ 
ous  publications  are  mostly  speculative,  and  assume  certain  simple  source  types 
in  the  ice  plate.  Also,  it  is  not  clearly  understood  yet  how  the  environmen¬ 
tal  stresses  are  released  as  a  major  source  of  induced  ice  motion.  The  possible 
candidates  are  the  three  dominant  types  of  crack  [18]  .  Since  the  field  obser¬ 
vation  of  such  cracks  does  not  seem  to  be  feasible,  modeling  of  radiation  from 
different  types  of  crack  is  comparative  to  infer  the  source  mechanism  from  the 
observed  signal  of  the  ice  event.  The  fracture  processes  are  expected  to  be  dif¬ 
ferent  depending  on  the  material  properties  as  well  as  the  stress  distribution,  of 
which  corresponding  radiation  patterns  can  be  characterized  by  a  set  of  source 
parameters,  such  as  fracture  type,  fault  orientation  and  dimension,  and  prop¬ 
agation  speed.  The  waveguide  nature  of  the  laterally  stratified  medium  affects 
the  radiated  field  as  well.  The  observed  signal,  therefore,  contains  informations 
characterizing  such  source  and  environmental  parameters.  The  result  of  this 
study  applied  to  source  inversion,  by  which  the  relationship  between  the  ob¬ 
served  acoustic  signal  and  the  physical  process  of  sound  generation,  will  lead  to 
better  understanding  of  source  mechanisms,  and,  consequently,  of  the  central 
Arctic  ambient  noise  as  an  aggregate  of  the  individual  event.  Another  expected 
application  can  be  found  in  the  development  of  remote  sensing  techniques  suit¬ 
able  for  the  study  of  ice  behaviour  and  properties. 
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Figure  1.1:  The  research  issue  applied  to  the  seismic  study  of  ice  and  the  ambient 
noise  generation  mechanism  due  to  ice  cracking  noise 
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1.2  Objective 


The  thrust  of  the  present  research  is,  therefore,  to  develop  a  model  of  the  elastic 
wave  field  in  stratified  environments  for  the  various  seismic  source  mechanisms. 
The  research  effort  in  this  thesis  is  focused  on  the  modeling  of  the  physical 
processes  of  radiation  from  given  types  of  source  and  propagation  in  a  range- 
independent  environment,  as  described  in  the  box  in  Fig  1.1,  i.e.  simulation  of 
the  radiation  patterns  and  time  series  for  various  given  types  of  sources.  The 
source  types  to  be  considered  are  explosive  sources,  point  forces,  shear  cracks, 
and  tensile  cracks.  Further  the  effect  of  transverse  isotropy  will  be  considered, 
which  is  important  in  sea  and  lake  ice,  and  finely  layered  sediments. 


1.3  Approach 

The  approach  taken  in  modeling  the  radiation  and  propagation  in  a  laterally 
stratified  medium  is  divided  into  two  part,  the  development  of  a  mathematical 
representation  of  seismic  sources  and  numerical  solution  technique  for  laterally 
stratified  elastic  media. 

The  representation  of  seismic  sources  utilizes  the  formulation  by  Keilis- 
Borok(1950),  reviewed  by  Sato[36],  which  introduces  the  displacement  and  force 
potentials  following  the  Love-Stoke  formalism[25] .  The  formulation  is  compact, 
and  can  be  easily  transformed  to  cylindrical  coordinates,  which  is  more  conve¬ 
nient  for  the  stratified  media  due  to  the  geometric  nature  of  horizontal  wave 
guide. 

The  methods  available  for  treating  the  acoustic  propagation  problems  can  be 
categorized  based  on  the  geometry  of  the  environmental  model  and  the  solution 
technique.  Confining  our  interest  to  the  seismic  and  acoustic  radiation  and  prop¬ 
agation  in  ocean  environments,  the  classification,  based  on  the  geometry  of  envi- 
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Figure  1.2:  The  solution  technique 

ronmental  model,  is  by  tradition  divided  into  two  categories,  range-independent 
and  range-dependent  media.  A  range-independent  medium  is  also  referred  to  a 
laterally  homogeneous,  laterally  stratified,  or  vertically  varying  medium. 

Since  range-dependent  model  allows  more  accurate  modeling,  a  great  amount 
of  effort  is  being  invested.  However,  the  solution  techniques  for  a  range-dependent 
medium  are  certainly  more  complicated  and  computationally  expensive  due 
to  the  more  complex  of  the  geometry.  The  available  methods  for  the  range- 
dependent  problems  have  certain  drawbacks.  For  example,  ray  theory  is  a  high 
frequency  approximation  with  singularities  such  as  caustics,  shadow  zone,  and 
head  waves,  which  need  special  treatment  [7].  Modal  approaches  such  as  adi¬ 
abatic  and  coupled  modes  for  the  range-dependent  environment  are  available, 
however  they  do  not  describe  the  near  field  properly.  Also,  the  available  codes 
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are  for  acoustic  media  only.  The  parabolic  equation  approach  is  applicable  to 
weakly  range-dependent  acoustic  media  only.  Numerically  intensive  methods, 
such  as  finite  difference  method  and  finite  element  method,  are  available,  but 
they  are  not  suitable  for  studying  the  long  range  propagation  for  computational 
reasons. 

Fortunately,  for  many  cases,  range-independent  models  can  be  applied  effi¬ 
ciently  for  many  ocean-sediment  environments  as  well  as  the  ice  covered  ocean 
found  in  the  central  Arctic,  except  when  it  is  necessary  to  consider  specific  fea¬ 
tures  such  scattering  due  to  ice  ridges.  This  class  of  models  can  be  modified 
to  treat  certain  small  imperfections  such  as  roughness.  Therefore,  the  solution 
techniques  available  for  range-independent  seem  to  be  more  suitable  for  our 
present  research  interest. 

The  most  complete  solution  technique  for  range-independent  environments 
is  the  integral  transform  methods.  By  means  of  temporal  and  spatial  integral 
transforms,  the  partial  differential  equation  of  motion  is  reduced  to  an  ordinary 
differential  equation  in  the  vertical  space  coordinate  z  (Fig  1.2.  The  inverse 
transforms  can  be  evaluated  analytically  by  means  of  asymptotic  approxima¬ 
tions,  such  as  steepest  descent,  or  stationary  phase.  However,  to  described  total 
fields,  the  inverse  transformation  must  be  evaluated  numerically.  The  transform 
method  can  be  classified  as  three  well  known  methods.  First,  matrix  method 
uses  the  time-frequency  Fourier  transform  in  order  to  formulate  the  equation  of 
motion  in  frequency  domain.  After  reducing  the  partial  differential  equation  to 
an  ordinary  differential  equation,  a  linear  system  of  equations  can  be  obtained 
for  each  discretized  horizontal  wave  numbers  to  evaluate  the  wave  number  in¬ 
verse  transform.  Second,  in  modal  method,  the  difference  to  matrix  approach 
is  that  the  integral  transform  is  replaced  modal  sum,  which  is  equivalent  to  ne¬ 
glecting  the  near  field  effects.  Also,  although  the  theory  is  being  developed  to 
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solve  for  elastic  media,  modal  method  is  mostly  limited  to  an  acoustic  medium. 
In  generalized  ray  theory,  the  time-frequency  inverse  transform  is  applied  first, 
then  the  inverse  wave  number  integration  is  evaluated  [6].  Although  ray  theory 
and  modal  method  can  be  extended  to  laterally  inhomogeneous  medium,  the 
matrix  approach  proves  to  be  the  most  widely  used  method  for  laterally  homo¬ 
geneous  medium.  Especially,  the  global  matrix  approach  applied  here  give  full 
wave  solution  for  all  layers  solving  the  azimuthal  Fourier  orders  simultaneously. 


The  solution  technique  in  a  laterally  stratified  medium  uses  cylindrical  co¬ 
ordinates.  The  advantage  of  using  cylindrical  coordinates  in  range-independent 
environment  is  that  only  one  integral  transform  in  variable  r  is  necessary  since 
the  azimuthal  angle  is  in  the  form  of  Fourier  summation,  while,  in  rectangular 
coordinates,  two  integral  transformations  in  x  and  y  are  involved  resulting  in 
the  numerical  inefficiency  when  evaluating  the  inverse  transforms.  In  cylindrical 
coordinates,  the  equation  of  motion  is  depth-separated  by  Fourier  transform  in 
the  azimuthal  angle  and  Hankel  transform  in  range.  The  remaining  ordinary 
differential  equation  in  the  depth  coordinate  z  with  proper  boundary  conditions 
for  horizontal  interfaces  yields  a  set  of  linear  system  of  equation  for  each  inter¬ 
face.  The  local  matrix  is  properly  combined  to  yield  the  global  matrix,  giving 
a  full  wave  solution  simultaneously  for  each  azimuthal  Fourier  orders,  and  for 
all  layers.  Once  the  depth-dependent  solutions  are  found  for  each  layers,  the 
inverse  Hankel  transform  is  performed  for  each  azimuthal  Fourier  order.  Now, 
the  frequency  domain  solution  is  found  by  summing  the  Fourier  orders.  Finally, 
the  time  domain  solution  is  found  by  frequency-time  inverse  Fourier  transform 
(Fig  1.2). 
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Figure  1.3:  Overview  of  the  thesis 

1.4  Overview  of  the  Thesis 

This  section  is  intended  as  an  introductory  overview  of  the  thesis.  The  topics 
treated  in  each  chapters  are  shown  in  Fig  1.3.  The  methodology  is  discussed  in 
somewhat  detail,  and  the  results  are  briefly  presented. 

In  Chapter  2,  the  compact  representation  of  seismic  sources  in  a  homoge¬ 
neous,  unbounded  medium  is  discussed.  First,  the  radiation  due  to  a  single  body 
force  with  arbitrary  direction  in  a  homogeneous  medium  is  derived.  Introduc¬ 
ing  the  Keilis-Borok’s  compact  formulation  of  arbitrary  source  time  function, 
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and  again  considering  harmonic  source  function,  the  formulation  in  frequency 
domain,  which  is  relevant  to  present  approach,  is  found.  Then,  the  modeling  of 
seismic  sources  in  a  homogeneous  medium  is  considered  through  the  represen¬ 
tation  theorem  and  Burridge  and  Knopoff  [5],  by  which  displacement  disconti¬ 
nuities  are  replaced  with  equivalent  body  forces.  Accordingly,  the  higher  order 
sources,  such  as  shear  crack  and  tensile  crack  sources,  are  represented  as  super¬ 
positions  of  force  couples  in  certain  directions.  In  order  to  construct  the  force 
couples  to  form  higher  order  sources,  a  unified  approach  to  represent  all  compact 
sources  by  force  potentials,  following  the  Love-Stoke  formalism,  is  discussed  and 
presented. 

In  Chapter  3,  in  order  to  incorporate  these  source  terms  into  a  propaga¬ 
tion  model  determining  the  full  solution  in  a  laterally  stratified  medium,  the 
depth-dependent  Green’s  functions  in  cylindrical  coordinates  are  derived  for 
each  type  of  source  in  Section  2.  The  homogeneous  solutions  in  a  laterally  strat¬ 
ified  medium  in  cylindrical  coordinates  are  assumed,  by  which  the  equation  of 
motion  is  separated  into  wave  equations  for  each  potentials  by  azimuthal  angle 
Fourier  transform  and  Hankel  transform  in  range.  The  remaining  ordinary  dif¬ 
ferential  equation  with  respect  to  z  with  a  proper  set  of  boundary  conditions  for 
horizontal  interfaces  yields  a  set  of  linear  system  of  equation  for  each  interfaces. 
The  local  matrices  are  properly  combined  to  yield  the  global  matrix,  giving  a  full 
wave  solution  simultaneously  for  all  layers.  Once  the  depth-dependent  solutions 
are  found  for  each  layers,  the  inverse  Hankel  transform  is  performed  for  each 
azimuthal  Fourier  orders.  The  frequency  domain  solution  is  found  by  summing 
the  Fourier  orders.  Time  domain  solution  is  found  by  frequency-time  inverse 
Fourier  transform.  In  this  research,  the  prototype  of  three  dimensional  version 
of  a  Fast  Field  Program  (FFP),  so  called  Seismo- Acoustic  Fast  field  Algorithm 
in  Range  Independent  environment  (SAFARI),  which  has  been  developed  by 
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Schmidt  and  Glattetre[38]  is  further  developed  to  take  seismic  source  terms, 
and  to  solve  the  propagation  in  a  stratified  medium.  This  particular  solution 
technique  is  chosen  due  to  the  efficient  simultaneous  treatment  of  the  several 
azimuthal  Fourier  orders  of  the  field. 

In  Chapter  4,  a  model  for  radiation  from  the  non-compact,  or  propagating 
crack  is  formulated,  first  for  an  unbounded  homogeneous  medium  and  then,  for 
a  stratified  medium.  The  field  caused  by  the  moving  crack  is  basically  found 
by  integrating  the  proper  source  Green’s  function  over  the  fault  surface,  where 
the  source  needs  to  be  placed  off  the  2-axis.  The  mathematical  treatment  of  the 
source  off  the  2-axis  involves  the  calculation  of  higher  order  Bessel  function  [38] , 
which  causes  the  numerical  inefficiency  and  convergence  problems.  Therefore, 
a  simple  transform  of  the  source  and  receiver  position  is  used  in  this  study  so 
that  the  source  is  always  on  the  2-axis.  Some  simple  examples  are  discussed  for 
which  the  analytical  solution  exists.  The  distinct  characteristics  of  the  radiation 
pattern  shows  the  directivity  pattern  dependent  on  the  ratio  of  wave  length 
to  crack  dimension.  When  this  ratio  is  large,  the  source  can  be  treated  as  a 
compact  source.  When  the  ratio  is  small,  the  field  is  highly  directional.  The 
synthetic  time  series  show  the  characteristics  of  the  convolved  source  signal 
with  the  dimension  of  the  crack  surface,  such  as  the  stopping  phase,  and  poor 
correlation  between  channels,  and  signal  duration  depending  on  the  observation 
positions. 

Next,  in  Chapter  5,  the  propagation  in  a  transversely  isotropic  medium, 
which  is  characterized  by  5  elastic  constants,  is  considered.  The  equation  of 
motion  is  no  longer  reduced  to  the  wave  equation.  However,  the  intermediate 
functions  can  be  defined  to  decouple  the  equation  of  motion  for  SH  and  SV  —  P 
waves  in  cylindrical  coordinates.  The  decoupled  equations  give  6  eigenvalues 
and  eigenvectors  for  up-going  and  down-going  components  of  SH  ,  SV ,  and  P 
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waves,  respectively.  Those  intermediate  functions  are  combined  to  formulate 
a  set  of  boundary  conditions  compatible  with  the  Global  Matrix  method.  An 
example  in  a  homogeneous  unbounded  medium  shows  the  expected  separation  of 
SH  and  SV  waves  when  propagating  in  horizontal  direction,  which  results  from 
the  different  shear  modulus  for  vertical  and  horizontal  directions.  The  slowness 
surfaces  are  given  and  discussed  for  the  example. 

Chapter  6  discusses  the  crack  radiation  from  the  3  modes  of  crack  in  a 
floating  ice  plate.  The  discussion  is  focused  on  the  characterization  of  each 
source  in  terms  of  radiation  and  temporal  characteristics.  The  radiation  from 
propagating  cracks  is  also  included  for  a  typical  set  of  source  parameters.  The 
variation  of  the  field  due  to  source  parameter  changes  are  briefly  discussed, 
including  the  effect  of  the  environmental  parameters  such  as  the  anisotropy  of 
the  sea  ice.  For  source  inversion  from  the  data,  more  importantly  to  test  the 
developed  model,  an  experiment  is  proposed.  The  objectives  of  this  experiment 
include  studying  the  relationship  between  the  observed  acoustic  signal  and  the 
physical  process  of  fracture  mechanism,  and  further  the  mechanical  behavior 
and  properties  of  ice  under  different  loading  conditions. 

Finally,  Chapter  7  summerizes  the  results  and  the  discussions  encountered 
through  the  research,  along  with  suggestions  for  further  studies. 
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Chapter  2 

Representation  of  Seismic 
Sources 


In  this  chapter,  the  representation  of  seismic  sources  in  an  unbounded  homoge¬ 
neous,  and  isotropic  medium  is  considered  in  terms  of  equivalent  body  forces. 

First,  the  representation  theorem  [l]  [4]  is  introduced  to  represent  the  fault 
motion  by  the  body  force  equivalents,  which  are  the  double  couples  [5],  In 
order  to  formulate  the  double  couples  in  a  compact  form  using  Keilis-Borok’s 
results,  the  radiation  from  a  single  force  is  discussed  in  detail.  The  force  and 
displacement  potentials  are  consistently  used  for  ease  and  well  defined  math¬ 
ematical  manipulation,  which  is  only  possible  for  an  isotropic  medium1.  The 
single  couples  are  shown  to  be  obtained  from  the  force  by  taking  proper  deriva¬ 
tives  depending  on  the  direction,  and  the  double  couples  are  found  from  the 
linear  combination  of  single  couples.  Finally,  the  body  force  equivalents  for 
seismic  sources  such  as  strike-slip,  dip-slip,  and  tensile  crack  with  certain  dip 
angle  are  summarized.  These  results  will  be,  in  Chapter  3,  transformed  into 

1The  formulation  of  homogeneous  solution  in  a  transversely  isotropic  medium  is  discussed  in 
Chapter  5 
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cylindrical  coordinates  in  order  to  allow  an  arbitrary  dip  angle,  yielding  a  depth 
dependent  source  functions. 


2.1  Representation  Theorem  and  Body  Force 
Equivalents  for  Seismic  Sources 

The  mathematical  description  of  seismic  sources  has  been  pursued  along  two 
different  lines [l]  :  First,  in  terms  of  body  forces  ;  second,  in  terms  of  disconti¬ 
nuities  in  displacement  or  strain  across  a  rupturing  fault  surface  or  across  the 
surface  of  a  volume  source.  When  the  representation  theorem  is  used,  it  can 
be  shown  that  the  discontinuity  of  displacement  and  stress  can  be  expressed  in 
terms  of  body  force  equivalent.  Consequently,  the  representation  theorem  was 
used  by  Burridge  and  Knopoff  [5]  to  show  that  the  body  force  equivalents  for 
fault  motion  are  double  couples.  The  details  will  be  given  in  the  followings. 

The  displacement  due  to  body  force  f,  displacement  discontinuity  [u],  and 
the  traction  [T(u,n)]  is,  from  the  representation  theorem,  expressed  as 

un{x,t)  =  f  drift  fp(ri,ri)Gnp{x,t-T;ri,0)dV(rj) 

I 

+  f  dr  f  f  {[ui(Z,T)}cijpqnjGnPiq(x.,t  -  r,  £,0)  (2-1) 

II 

-[Tp(u(e,r),n)]Gnp(x,f-r;e,0)}dE(e)  . 

v . — . . . .  V - - 

III 

The  first  term  represents  the  body  force  contribution,  while  II  and  III  are  due 
to  the  displacement  and  traction  discontinuity  contributions,  respectively.  Now, 
neglecting  the  presence  of  body  force  and  traction,  only  the  second  term  remains. 
And,  using  the  relation 

-^-Gnp{x,t  -  r;  £,0)  =  -  J j ~  0GnP{x,t  -  T-Tj,0)dV(r})  ,  (2.2) 
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Figure  2.1:  Fault  motion. 

the  second  term  in  Eq  2.1  due  to  displacement  discontinuity  reduces  to 

un{x,t)  =  Jjlr  J I fln](r},T)Gnp(x,t  -  T]rj,0)dV(ri)  ,  (2.3) 

where  the  equivalent  body  force  for  the  displacement  discontinuity  [u]  is 

=  ~  f  Jjui{C,r)}cijpqn~6{v  -  £)dE  .  (2.4) 

Note  that  the  surface  integral  for  displacement  discontinuity  in  the  second  term 
in  Eq  2.1  has  been  expressed  in  terms  of  body  force  equivalent  in  Eq  2.4,  of 
which  form  is  found  in  the  first  term  I  in  Eq  2.1. 

As  an  example  of  a  buried  fault  as  shown  in  Fig  2.1,  the  equivalent  body 
forces  can  be  derived  from  Eq  2.4, 

Mv,r)  =  -M06{nx)S{rh)4-6MXo{T)  , 

/2(»?,r)  =  0,  (2.5) 

fs[r],T )  =  -M0-^%i)%2)<$(»?3)Aro(r)  , 
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▲ 


(a) 


(b) 


Figure  2.2:  Body  force  equivalents  for  derivatives  of  Green’s  function. 

where  the  seismic  moment  M0  in  an  isotropic  medium  is  defined  as 

Mq=huA=h  x  average  slip  x  fault  area  .  (2.6) 

The  corresponding  displacement,  from  Eq  2.3  using  Eq  2.2,  is  expressed  as 

«„(*,()  =  /_>//emW  dE.  (2.7) 

In  fact,  the  first  and  second  terms  in  the  curly  bracket  in  Eq  2.7  represent  single 
couples  in  Fig  2.2.  Since  the  two  forces  are  combined  in  such  a  way  that  the 
magnitude  of  moment  remains  constant  while  e  goes  to  zero,  i.e.  limt_l0  eX0  = 
M0,  the  single  couples  appear  to  be  derivatives  of  a  single  body  force. 

In  the  following  sections,  the  more  compact  form  of  seismic  source  representa¬ 
tion  using  displacement  and  force  potentials  is  derived  based  on  the  Love-Stokes 
formalism,  following  Sato  [36], 
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2.2  Radiation  from  Body  Forces  in  a  Homoge¬ 
neous  Isotropic  Medium 


As  a  basic  ingredient  in  the  representation  of  the  various  sources,  the  radiation 
field  due  to  a  single  body  force  is  derived  in  detail. 

The  displacement  equation  of  motion  for  a  homogeneous  isotropic  medium 
is 


MV’u  +  (A  +  ,<) V  (V  •  u)  +  pF  =  p-gjj  , 

(2.8) 

where  n  and  A  are  Lam4  constants  and  u  is  displacement,  p  is 

density,  and  F 

is  a  body  force.  Or,  equivalently 

(A  +  2p)  V  (V  •  u)  —  x  V  x  u  +  pF  =  p . 

(2.9) 

Taking  the  divergence  and  curl  of  Eq  2.9,  we  obtain 

cc2V2V  •  u  +  V  ■  F  —  gt,  , 

(2.10) 

and 

2_2r,  _  _  d2  V  x  u 

c,2V2V  x  u  +  V  x  F  =  at. 

(2.11) 

with 

(X  +  2fi\1/2  ( 

Ce  =  (  /-  J  ’c*  =  ui  ■ 

(2.12) 

being  the  compressional  and  shear  velocities,  respectively. 

Introducing  body  force  potentials  and  displacement  potentials, 


F  =  (X,  Y,  Z)  =  V$  +  Vx(i,  M,  N) ,  (2.13) 

u  =  (u,v,w)  =  V<f>  +  V  x  (F,  G,  H) ,  (2.14) 

Eqs  2.10  and  2.11  become 


c02VV  +  $ 


d2<t> 
dt 2 
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(2.15) 


c,2V2F  +  L  = 
c2V2G  +  M  = 
c2V2H  +  N  = 


d2F 
~di*~ 
d2G 
dt 2 
d2H 
dt 2 


Now,  the  particular  solution  of  Eq  2.15  can  be  expressed  in  the  forms. 


d>  = 

t  — —  j  dx^dx^dx's 

F  = 

5J.///M' 

f  — —'j  dx^dx^dxg 

(2.16) 

G  = 

4 

t - ^  dx^dx'zdx'z 

H  = 

4.W//M 

t  — —  J  dx^dx^dx^ 

ct/ 

where  $,L,M  and  N  are  found  from 


$ 

L 

M 

N 


dx^dx'zdx'z 


dx\dx\dx's 


dxxdx2dxz 


dx^dx'zdx's 


(2.17) 


where  t  is  the  time  and  r  =  \J (zi  —  x'y)1  +  (x2  —  x‘2)2  +  (X3  —  x'3)2.  The  proof  of 
Eq  2.16  for  displacement  potentials  is  drawn  basically  from  the  Possion  equation 
solution  in  an  unbounded  medium  [25,  pp.304]  [l,  pp.64-67]  [35,  §276].  Eq  2.17 
for  force  potentials  is  also  proved  utilizing  the  solution  for  vector  Poisson  equa¬ 
tion  2  [25,  pp.184]  [1,  pp.69].  In  the  above  expressions,  the  unit  of  vector  quantity 

2  Since  the  proof  is  elegant  and  educational,  it  is  given.  Force  is  defined  through  force  poten¬ 
tials,  F  =  (X,  Y,  Z)  =  V4>+ V x  (L,  M,  N).  An  vector  Poisson  equation  for  W  can  be  constructed 
with  F  as  source  term,  F  =  V2W,  which  can  be  rewritten  as  F  =  V(V  •  W)  —  V  x  (V  x  W). 
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X'  is  acceleration,  i.e.  m/sec2,  and  the  scalar  quantity  X0(t  —  r/ce)  is  source 
time  function  and  has  the  unit  of  force.  Finally,  substitution  of  Eq  2.17  into 
Eq  2.16  yields  the  displacement  potentials.  The  volume  integral  for  displace¬ 
ment  potentials  can  be  evaluated  for  an  arbitrary  source  function  by  conversion 
into  a  surface  integral  ;  details  are  given  by  Aki[l,  pp.  71]. 

The  general  formulation  of  particular  solution  to  the  wave  equation  with 
a  body  force  source  in  an  unbounded  medium  has  been  derived,  where  the 
source  function  can  be  completely  arbitrary.  An  example  for  a  horizontal  force 
with  haxmonic  time  dependence  will  be  derived  following  the  steps  presented 
previously.  Later  in  Chapter  2.5,  when  the  body  force  equivalents  for  seismic 
sources  are  formulated  in  the  form  of  force  potentials,  the  results  will  be  used 
to  find  displacement  field  due  to  corresponding  force  potentials.  The  usage  of 
harmonic  source  function  has  been  chosen,  because  the  time  domain  solution  for 
an  arbitrary  source  function  can  be  found  by  time-frequency  Fourier  synthesis. 

Example  :  Radiation  from  a  single  force  in  x-direction 

For  the  case  of  horizontal  force  in  xi-direction  with  harmonic  time  dependence 
i.e.  F  =  X06(x)ei,  Eq  2.17  is  simplified  as  follows, 

$  =  ~4 nJIl{X'^)dx'ldX2dx'3 

L  =  0  (2.19) 

m  = 

It  is  noted  that  the  force  potentials  are  now  expressed  in  terms  of  W  as,  $  —  V  •  W,  and 
(L,  M,  N)  —  V  X  W.  The  solution  W  for  the  vector  Poisson  equation  is 

w=-///  4^er(f)  (2i8) 

When  the  above  equation  2.18  is  substituted  into  the  expressions  for  force  potentials,  $  =  V  -  W 
and  (L,  M,  N)  =  V  x  W,  Eq  2.17  is  proved. 
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N  -  4tt  ///  (X>  dx2  )  dx'idx*dx'* 

The  volume  integrals  in  Eq  2.19  can  be  found  using  relation 

P  J J J  X' dx[dx'2dx'3  =  X0  ( t ) ,  (2.20) 

where  ei  is  an  unit  vector  in  £i-direction,  and  X  —  XqS(x'),  and  the  scalar 
quantity  Xo{t )  is  a  source  time  function.  Then,  the  force  potential  reduces  to 


$  = 


X0  dr 


-l 


4 irp  dx\ 


L= 0 ,  M = 


X0  dr-1 


M  = 


Xo  dr 


-l 


- 5—  ,  -  -  •  (2.21) 

47rp  ox 3  47rp  dx2 

Substituting  Eq  2.49  into  Eq  2.16,  the  displacement  potentials  are  found  to  be 

=  “ (idjV  ///  (*  “  7)  %^dx'idx2dx> 

F(r,t)  =  0  (2.22) 

G(r,()  =  (4ttc,)v  Iff  7ro  (‘ _  t)  -a^d<dx'iix* 

The  volume  integral  for  displacement  potentials  can  be  evaluated  for  an  arbitrary 
source  function  by  conversion  into  a  surface  integral  ; 

i  d  1  f"  TX^t  _  Tjdr 

4irp  oil  r  Jo 


<f>  = 


F  =  0 
G  = 


i  d  i  rh 


(2.23) 


47 Tp  dx  3 

H  =  1  3  1  /■* 


—  J°’  rXo(t  —  r)dr 


TX0(t-r)dT 

4np  ox 2  r  Jo 
details  are  given  by  Aki[l,  pp.  71]. 

Evaluation  of  the  volume  integral  in  Eq  2.22  given  by  Keilis-Borok  (1950), 
reviewed  in  Sato  [36],  gives  a  compact  result,  replacing  the  volume  integrations, 
by  introducing  the  intermediate  functions  defined  as 


,  1  r .  ,  1  .  r  . 

4>  o  =  ~F(t - )  ,  tfo  =  -F(t - ) 

r  ce  r  c. 


(2.24) 
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(2.25) 


F(t)  =  f*ds'  ['  X0  {s)d* 

Jo  Jo 

where  Xq (t)  is  a  source  function.  For  a  harmonic  time  dependence  of  source 
function,  i.e.  X0(i)  =  X0e,ult,  the  Keilis-Borok’s  intermediate  functions  0O  and 
0o  are 

M*)  =  -~e^‘-r/ec'>  ,  0o  (0  =  (2.26) 

urr  u2r 

Using  this  notation,  the  volume  integrations  for  displacement  potentials  in  Eq  2.22 
reduce  to 

^  =  1  d<f>Q  jr  -q  G=. _ 1  d^o  s  _  1  dTpo 

4i rp  dxi  ’  ’  47rp  dxs  ’  47rp  dx2 

It  is  pointed  out  that  Eq  2.27  is  equivalent  to  the  result  by  surface  integration  in 
Eq  2.23.  Using  the  relation  given  in  Eq  2.14  for  the  displacement  components, 
we  obtain 


u  =  V0  +  V  x  (F,  G,  H ) 

=  V0  +  V  x  V  x  (0i,02, 03)  (2.28) 

where 

0  =  — — 00?  01  =  —  7— 0o>  02  =  0,  03  =  0  (2.29) 

47rp  47rp 

Similar  expression  may  be  found  easily  when  the  forces  are  applied  in  x2  and 
z3-directions.  The  radiation  field  due  to  a  force  in  an  arbitrary  direction  also 
can  be  found  by  superposition  of  the  decomposed  force  components  in  the  xi,x2, 
and  x3-directions. 

The  above  results  are  remarkably  simple,  and  will  be  used  subsequently  to 
formulate  the  field  expressions  for  higher  order  sources. 
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Figure  2.3:  Decomposition  of  single  couple  into  two  body  forces. 

2.3  Single  Couples 

Potentials  for  single  couples  can  be  found  by  combining  the  single  forces  previ¬ 
ously  formulated.  Instead  of  the  force  source  function  X0(t),  the  seismic  moment 
M0{t)  is  used  hereafter  for  representation  of  higher  order  sources,  so  that  the 
definitions  of  <f>o  and  ipo  are 

<£0  =  iF(i--)  ,  V’o  =  ~F{t  -  — )  (2.30) 

v  cc  r  Cj 

F{t)  =  r  ds'  r  M0(s)ds  (2.31) 

Jo  Jo 

where  Mo(t)  is  a  moment  source  function. 

The  first  single  couple  in  the  Fig  2.5,  for  example,  can  be  generated  by 
combining  two  xi-direction  single  body  forces,  of  which  formulations  are  given 
previously.  The  displacement  potentials  for  the  single  couple  in  the  Fig  2.3  are 
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easily  found  to  be 


^  =  1  ^Vo 

47 xp  dxzdxi 

F  =  0 


Q  _ _ 1  d*rl>o 

47r p  dxl 

L  =  i  ay0 

47rp  dx2dx3 

The  displacement  can  now  be  written  in  compact  form  as 


(2.32) 


u  =  V^  +  VxVx  (0x,O,O)  , 


(2.33) 


where  the  potentials  are 


.  _ _ 1  dVo  _  1  dt[>0 

4xp  dxzdxi  ’  1  47rp  dx3 


(2.34) 


For  the  second  single  couple  in  Fig  2.5  with  forces  parallel  to  the  x3-axis,  the 
potentials  are 


4>  - 


1  d2<f>0 


47 xp  dxidx3 

F  _ _ 1  d2rp0 

4 ixp  dx2dx\ 
G  _  1  d2rp o 

47rp  dx\ 

L  =  0  . 


G  = 


(2.35) 


Again  using  a  compact  form,  the  displacement  is 


u  =  V<f>  +  V  x  V  x  (0, 0,  tps) 


(2.36) 


where  the  potentials  are 


.  _ 1  d2<j> o  _  1  drpp 

4'Kpdx\dxz  ’  3  47rp  5xx 


(2.37) 
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These  single  couples  are  exactly  same  as  the  example  (Fig  2.2)  discussed  in 
Chapter  2.1  with  different  notations.  The  single  couples  in  Eq  2.34  and  Eq  2.37 
correspond  to  the  first  and  second  terms  in  the  curly  bracket  in  Eq  2.7. 

Potentials  for  the  nine  single  couples  (Fig  2.4)  that  are  required  to  obtain 
equivalent  body  forces  for  a  generally  oriented  displacement  discontinuity  are 
summarized  below. 


(1,1) 

(1,2) 

(1.3) 
(2,1) 
(2,2) 

(2.3) 

(3.1) 

(3.2) 

(3.3) 


<f> 

<P 

<t> 

<f> 

4> 

4> 

<f> 

4> 

<f> 


1  d24>o  ,  \8ip0 

=  T— 7T~,  ^2  =  0,  ip  %  —  0 


4t xp  8x\  ’ 

i  ay0 

4 irp  8x38x1  ’ 
1  d2<p0 

4irp  dxsdxi  ’ 

1  dVo 


V’i  = 


47 Tp  dll 
1  dipo 


4np  dx2  ’ 


^2  =  0,  ip  3  —  0 


47rp  dxidx2  ’ 
1  d2<po 


i  1  dtpo  , 

Wl  =  -A - ~ - ,  V>2  =  0,  tp3  —  0 

47rp  ax3 

1  drpo 

i>i  =  0,  xp2  =  -—^,  ips  =  0 


4irp  dx\  ’ 
1  d24>0 


rpx  =  0,  xpt 


47 xp  dxi  ’ 

1  dipo 


4tx p '  dx2  ’ 


^3  =  0 


47rp  8x38x3  ’ 

1  ay0 

4717?  8x18x3  ’ 

1  ay0 

4lXp  8X38X3  ’ 

1  5Vo 


=0,  Ip  2  =  -7—^—,  ^3  =  0 


47rp  8x3  ’ 

V'i  =  0,  V>2  =  0,  ^3  = 

Ipi  =  0,  ^2  =  0,  rp3  = 


1  d^p 
47rp  8x  1 
1  5^0 


47Tp  $£3  ’ 


Vh  =0,  =  0,  xps  = 


4ixp  8x2 

1  dtpo 


4i xp  8x3 


(2.38) 

(2.39) 

(2.40) 

(2.41) 

(2.42) 

(2.43) 

(2.44) 

(2.45) 

(2.46) 


Using  the  single  couples,  the  expression  for  a  moment  are  found.  For  exam¬ 
ple,  the  x— component  for  a  moment  is  a  superposition  of  (3,2)  and  negative 
of  (2,3).  As  expected,  the  compressional  part  vanishes,  while  only  the  shear 
potentials  remain. 


41 


42 


Figure  2.6:  Strike-slip  with  dip  angle  6  =  0  (a)  fault  surface  and  slip  motion  (b) 
equivalent  body  forces. 


2.5  Summary  of  Equivalent  Body  Forces  for 
Seismic  Sources 


Strike-slip 


The  double  couple,  which  is  an  equivalent  body  force  for  strike-slip  with  dip 
angle  6  =  0  (Fig  2.6),  is  simply  a  superposition  of  two  single  couples  (1,3)  and 
(3,1),  giving  the  following  displacement  potentials, 


^  _ _ 1  d24>0 

2 wp  dxi  dx3 


= 


47 rp  dxs 


^2  =  0,  xJj3  = 


47rp  dx\ 


(2.49) 


Dip-slip 

The  displacement  potentials  for  dip-slip  (Fig  2.7)  with  dip  angle  6  =  0,  is  a 
superposition  of  two  negative  single  couples  (2,3)  and  (3,2)  in  Fig  2.4,  where 
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Figure  2.7:  Dip-slip  with  dip  angle  6  =  0  (a)  fault  surface  and  slip  motion  (b) 
equivalent  body  forces. 


the  potentials  axe 


27 xp  dxzdxz  ' 


n  ,  1  1  dtpo 

0  ,  ^2  =  —T—^—  >  V>3  =  •  2.50 

47 rp  ox 3  47rp  c/X2 


The  representation  of  shear  seismic  sources  in  rectangular  coordinate  system 
has  been  reviewed  for  certain  direction  of  fault  surfaces  using  Keilis-Borok’s 
compact  results.  The  more  complete  formulations,  using  the  nine  single  couples, 
in  rectangular  coordinates  for  general  shear  fault  surfaces  with  fault  orientation 
parameters(Fig  3.2)  such  as  dip  angle (<5),  rake  angle(A),  and  strike^),  are 
found  in  Aki  [1,  pp  117-118]  .  The  corresponding  representation  in  cylindrical 
coordinates  will  be  given  in  Chapter  3.2. 


Tensile  crack 

The  equivalent  body  force  for  a  tensile  crack  can  be  obtained  by  superposing 
three  single  couples  without  moment,  i.e.  (1,1),  (2, 2),  and  (3,3)  in  Fig  2.4,  but 
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Figure  2.8:  Tensile  crack  with  dip  angle  6  =  0  (a)  fault  surface  and  motion  (b) 
equivalent  body  forces. 


with  different  magnitude  of  seismic  moment.  Considering  the  fault  surface  with 
dip  angle  6  =  0  (Fig  2.8),  the  direction  of  fault  motion  is  in  ^-direction,  so  that 
the  magnitude  of  the  single  couple  (3,3)  will  be  ^^Mo,  while  the  other  two 
components  are  Mq.  The  potentials  for  tensile  crack  shown  in  Fig  2.8  are  then 


<t>  =  ~ 


i  av o 


1  d2&)| 


i  aVo 


*  i  a _ r\ 


4ttp  dx\  |M=Mo  4t rp  dx\  \M=Mo  47 rp  dx\  |m=a^Mo 
1  dtpo  I  , 


V»2  = 


= 


4np  dXl\MzzMo 

1  dtpo 

4-np  dx2  M=Mg 

1  dtpp 

dx 3 


(2.51) 


The  potentials  for  explosive  source  can  be  found  by  taking  same  magnitude 
in  3  single  couples  without  moment,  yielding  only  compressional  potentials. 
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Chapter  3 

Field  Representation  in  a 
Laterally  Stratified  Medium 


Previously,  in  Chapter  2,  the  displacement  field  in  a  homogeneous  isotropic 
medium  excited  by  various  kind  of  sources,  such  as  forces,  single  couples,  mo¬ 
ments,  double  couples  which  corresponds  to  shear  fault  motion,  and  body  force 
equivalents  to  tensile  crack  opening,  has  been  discussed  in  Cartesian  coordi¬ 
nates.  In  this  Chapter,  the  field  representation  in  a  laterally  stratified  medium 
is  considered  based  on  global  matrix  method.  First,  the  homogeneous  solution  in 
cylindrical  coordinates  is  presented  in  Section  3.1,  since,  in  a  range-independent 
medium,  two  integral  transforms  in  x  and  y  in  rectangular  coordinates  are  re¬ 
duced  to  a  single  integral  transform  in  r  in  cylindrical  coordinates  Then,  in 
Section  3.2,  the  displacement  potentials  in  the  previous  section  will  be  trans¬ 
formed  to  represent  the  dislocation  sources  with  arbitrary  dip  angle  <5  =  0. 
These  potentials  are,  again,  transformed  into  cylindrical  coordinates.  In  cylin¬ 
drical  coordinates,  strike  angle  (<j>s  in  Fig  3.2)  can  be  accounted  for  by  rotation  in 
azimuthal  angle.  Since  the  homogeneous  and  source  terms  are  in  the  same  form 
of  representation  via  Fourier  decomposition  in  azimuthal  angle  0,  Hankel  inte- 
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gral  transformation  in  r,  and  finally  the  depth-dependent  exponential  form  with 
unknown  coefficients,  a  set  of  linear  equations  can  be  established  by  matching 
boundary  condition  at  interfaces,  for  each  discretized  horizontal  wave  number  of 
Hankel  integral  representation  and,  for  each  Fourier  order  of  0,  which  is  treated 
in  Section  3.3.  In  Section  3.4,  the  solution  field  parameters  of  interest  such 
as  stresses  and  displacements  are  found  by  inverse  Hankel  transform  and  then, 
adding  all  the  Fourier  components.  Also,  the  numerical  aspects,  concerning  the 
inverse  Hankel  transform  are  discussed.  Finally,  some  numerical  examples  are 
given  in  Section  3.5. 


3.1  Homogeneous  Solution  to  Wave  Equation 
in  Cylindrical  Coordinates 


The  homogeneous  solutions  to  wave  equations  in  cylindrical  coordinates  are  rep¬ 
resented  in  terms  of  three  scalar  potentials,  which  are  related  to  the  displacement 
in  the  following  way  [38], 


u  =  V<£  +  V  x  V  x  (0,0,  A)  +  V  x  (0,0,  ip)  (3.1) 

The  scalar  potentials  satisfy  the  following  homogeneous  wave  equations, 

(V2  +  h2)<f>  =  0  (3.2) 

(V2  +  A2)(A,^)=0  (3.3) 

The  potentials  are  expanded  in  Fourier  series  in  the  azimuthal  angle  0,  as 

cos  m0 


4>{r,0,z)  =  <f>m{r,z) 

m= 0 


A  (r,0,z)  =  A  m(r,z) 

m=  0 


sin  mO 

cos  m0 
sin  m0 


(3.4) 
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oo 

^M.*)  =  Y  Vm{r,z) 

m= 0 

By  insertion  of  Eq  3.4  to  Eq  3.1,  it  can  be  shown  that  the  displacement  compo¬ 
nents  are  similarly  expanded  as 


w(r,0,z)  =  Y  wm(r,z) 

m=  0 


cos  m9 
sin  mO 


sin  mO 
—  cos  mO 


u{r,e,z )  =  Y,  «m(r>2) 
m=o  • 


cosm# 
sin  mQ 


(3.5) 


v(r,6,z)  =  '£vm(r,z) 

m= 0 

Substitution  of  Eq  3.4  into  the  wave  equations  3.2  and  3.3  results  in  ordinary 
differential  equations  for  the  expansion  coefficients  in  Eq  3.4  after  Hankel  trans¬ 
form  in  range  [38].  The  ordinary  differential  equation  with  respect  to  z  is  solved 
to  give  the  exponential  form  in  the  square  bracket  in  the  following  equation  3.6, 
which  represent  the  up-going  and  down-going  waves,  with  the  following  integral 
representations  for  m-th  order  potentials. 

<f>m(r,z)  =  J~  [ar(s)e-*“W  +  a?(s)e‘“M]  sJm{rs)ds 
A m(r,z)  =  J”  [6r(s)e-^  +  b?(s)e*^]  Jm(rs)ds 
ipm{r, z)  =  J™  [c? (s)e-^W  +  c«(a)e^W]  sJm{rs)ds 

In  order  to  match  the  boundary  conditions  at  the  interface,  it  is  necessary  to 
find  the  field  parameters  such  as  displacements  and  stresses.  The  constraints  at 


(3.6) 


sin  mQ 
—  cos  mQ 
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the  horizontal  interfaces  between  2  solid  media  in  welded  contact  involve  : 


w  :  z  -direction  displacement 
u  :  r  -direction  displacement 
v  :  0  -direction  displacement 
ozz  :  normal  stress  in  z  -direction 
orz  :  shear  stress  in  r  -direction 
Gfz  :  shear  stress  in  0  -direction. 


When  one  or  two  of  the  media  joining  at  an  interface  are  liquid  or  vacuum, 
some  of  the  constraints  will  be  no  longer  effective.  These  boundary  conditions 
are  used  to  determine  the  unknown  constants  which  are  functions  of  horizontal 
wave  number  s ,  such  as  aj(s),  a2(s),  &i(s),  &2(s)>  ci(5)  and  c2 Ca¬ 
using  Eq  3.6  and  the  relation  between  the  displacement  and  displacement 
potentials,  Eq  3.1,  the  field  parameters  given  above  can  be  shown  to  be  [38] 


roc 

wm(r,z)  = 

Jo 

um(r,z)±vm{r,z)  =  [ °° 

Jo 


-af(s)a(s)e-iaW  +  a™(s)a(s)e*“M 

+  b?(s)sezpW 


sJm(rs)ds 


T  a?{s)setaW 

±b?(s)0(s)e~zf)W  T  b?{s)0(s)ezl,W 

+c?{s)se-zl3M  +  c?(s)sezl3W 


sJm±i{rs)ds 


aTz(ri z)  =  A  VV*(r,z)  +  2»—wm{r,z) 


=  »./■ 


af(s)(2s2  -  A:2)e"*“W  +  a^(s)(2s2  -  k2)ezaW 
-b^{s)2s0(s)e-z^  +  b^(s)2s0(s)ez^ 


sJm(rs)ds 


(3.7) 


a™{r,z)  ±oZ{r,z)  =  n 


Yz  [um{r,z)  ±  vm{r,z)}  +  =F  \  wm(r,z ) 
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5^m±l(5r)d5 


±a™(s)2sot($)e  ^  a^s^safs)^01^ 

Tb?{s)(2s*  -  k2)e~zf)M  =F  6jT(s)(2s2  -  k2)ez^ 
-c?{s)sP{s)e-zM  +  c™(s)s/?(s)e^W 


In  the  above  equation  3.7,  it  is  noted  that  the  shear  displacements  and  stresses 
are  added  and  subtracted  giving  the  Hankel  transform  for  a  single  order  of  Bessel 
function,  i.e  Jm+1  for  addition  and  Jm_i  for  subtraction,  respectively.  Otherwise, 
the  shear  displacements  and  stresses  in  certain  directions  will  involve  Bessel  func¬ 
tions  of  two  different  orders,  which  greatly  complicates  analysis  using  a  linear 
system  of  equations.  This  is  an  essential  manipulation  in  order  to  be  able  to  ap¬ 
ply  global  matrix  approach.  The  homogeneous  solutions  to  the  wave  equations, 
and  corresponding  displacements  and  stresses  have  been  given.  The  coefficients 
are  later  found  by  matching  the  boundary  conditions  for  a  given  excitation.  In 
the  following  section,  the  source  potentials  are  transformed  into  the  cylindrical 
coordinate  system  to  have  same  form  of  representation  as  homogeneous  solutions 
found  in  Eq  3.6. 


3.2  Seismic  Source  Green’s  Function  in  Cylin¬ 
drical  Coordinates 

The  displacement  potentials  for  various  sources  given  previously  needs  to  be  for¬ 
mulated  in  cylindrical  coordinates  with  the  same  form  of  Fourier  decomposition 
in  ^-direction  as  applied  to  the  homogeneous  solutions  in  order  to  set  up  a  global 
linear  system  of  equations  expressing  the  boundary  conditions.  For  strike-slip, 
the  derivations  are  given  rather  in  detail,  while  only  the  results  are  given  for  the 
dip-slip  and  tensile  crack. 
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(a) 


(b) 


(c) 


Figure  3.1:  Fault  motion  with  dip  angle  8 ,  (a)  Tensile  crack,  (b)  dip-slip,  (c) 
strike-slip. 

3.2.1  strike-slip 


The  displacements  in  terms  of  displacement  potential  given  in  Eq  2.49  will  be 
first  transformed  into  a  strike-slip  with  dip  angle  8. 

Denoting  the  original  coordinates  system  (ii,  z2,  £3)  and  using  a  following 
transformation  (Fig  3.1),  the  strike-slip  has  dip  angle  8  in  the  coordinate  system 
{x,y,z), 


<f>  = 

0z  = 

0y  = 

4>z  - 


1  d 
2np  dx 
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smi- — bcosd— 
dy  dZj 


4np 


.  J  ,  xd' 

sin  b— — b  cos  0  — — 
dy  dz 


4zp 


sin  8 
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4'Kp 


cos  8 


dx 
dx  ’ 


00 


(3.8) 
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Again,  the  potentials  in  the  rectangular  coordinate  system  ( x ,  y,  z)  needs  to  be 
transformed  into  cylindrical  coordinate  system  (r,  0,  z).  Note  that  the  relations 
between  the  displacements  in  two  coordinate  systems,  for  example,  are 


ux  —  ur  cos  0  —  ug  sin  0 

uy  =  ur  sin  0  +  u$  cos  0  (3.9) 

uz  —  uz 


The  partial  differential  operators  can  be  found  using  the  chain  rule, 

d  sin  0  d  \ 


d_ 

dx 


=  ^COS0 


dr 


d  (  .  nd 

Ty=[SmeYr  + 


r  80t 
cos0  d 


80 


(3.10) 


d_  _  d_ 
dz  dz 

Using  the  relations  in  Eq  3.9  and  3.10,  the  potentials  are  found  as  follows.  The 
displacement  expression  equivalent  to  Eq  3.8  in  cylindrical  coordinates  is 


u  =  V^+VxVx  , 


(3.11) 


where 


4>  = 
Vv  = 
0*  = 


47rp 

1 

4  irp 


.  p  ( d2  1  8  \  82 

sin  o  sin  20  -—-r - —I  —  2  cos  8  cos  0---- 

\  8rl  r  dr  I  drdz 


<t>  o 


8  8 

—  sin  8  sin  20—  +  cos  8  cos  0— 
dr  dz 
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(3.12) 
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sin  6  cos  20  +  2  cos  8  sin  0 


dz 
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i  1  a 

ipz  =  - —  COS  0  cos  0- 


47 rp  dr 

The  potentials,  <f>0  and  t/)0,  can  now  be  expressed  as 


4>  o 


_  M°  1  -iutt-R/e.) 

w2  R 
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(3.13) 


=  -^-eiwt  f°°  J0{sr)e-alz-z-\-ds 
u2  Jo  a 


rpo  =  -^-etut  f  J0(sr)e  ^  z’^ds  , 
or  Jo  0 


using  the  Sommerfeld-Weyl  integral  [45,  pp.  13]  , 


a-ihR 


R 


roo  a 

=  /  J0(s»')c-a|2:l— ds  , 

Jo  a 


(3.14) 


where  a  =  ( h 2  —  s2)1/2fors2  >  Re(h2)  ,  and  j(h?  —  s2)1/2fors2  <  Re(h2),  J0 
is  zeroth  order  Bessel  function,  h  is  medium  wave  number,  R  is  y/r2  +  z2  , 
and  s  is  horizontal  wave  number.  It  is  noted  that  z  in  Eq  3.13  is  replaced  by 
z  —  zt  to  account  far  the  source  depth  zt,  accordingly  the  definition  of  R  is 
R  =  yjr*  +  (z  -  Za) 2. 

Finally,  substitution  of  the  above  Eq  3.13  into  Eq  3.12  yields  the  potentials 
of  strike-slip  with  an  arbitrary  dip  angle  6  of  the  same  integral  form  as  the 
homogeneous  solution  (note  that  f  is  defined  as  sign(z  —  zt)  in  the  following 
expressions)  : 
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Unlike  the  homogeneous  solution,  the  potentials  4>,fpr,4’e,fpz  rather  than  4>,k,4> 
are  used  in  particular  solution,  since  the  particular  solutions  are  found  by  trans¬ 
formation  of  the  rectangular  coordinate  formulation  to  give,  naturally,  the  cou¬ 
pled  potential  representation,  while  the  uncoupled  homogeneous  solutions  are 
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found  from  the  separated,  uncoupled  wave  equations.  The  expressions  in  the 
same  form  as  the  homogeneous  solution,  0,  A,  and  0,  can  be  shown  to  be 
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47rpo;2 

For  the  purpose  of  formulating  the  boundary  conditions  at  the  interface,  the 
displacements  and  stresses  can  be  calculated  from  either  of  potentials  in  Eq  3.15 
and  3.16. 


3.2.2  Dip-slip 

The  potentials  for  dip-slip  are  similarly, 
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(3.17) 


Referring  to  Fig  3.2  ,  strike-slip  corresponds  to  the  fault  motion  with  the  rake 
angle  A  =  0°  ,  and  dip-slip  to  A  =  90°.  Therefore,  the  fault  motion  with  a  rake 
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Figure  3.2:  Definition  of  fault  orientation  parameters. 

angle  A  is  a  superposition  of  the  potentials  for  strike-slip  and  dip-slip  multiplied 
with  cos  A  and  sin  A  ,  respectively.  The  strike  angle  can  be  incorporated  easily 
in  cylindrical  coordinates  by  replacing  the  azimuthal  angle  0  by  8  —  <ps  . 

It  has  been  shown  that  the  general  shear  fault  motion  in  cylindrical  coordi¬ 
nates  can  be  formulated  by  simple  coordinate  transformation.  Same  procedure 
will  be  applied  to  tensile  crack. 

3.2.3  Tensile  Crack 

The  tensile  crack  in  Fig  3.1(a)  is  a  superposition  of  three  single  couples.  Since 
the  cylindrical  coordinates  are  used,  the  potential  of  a  single  couple  in  Fig  3.3(a) 
can  be  used  to  express  the  remaining  two  single  couples  by  coordinate  transfor- 
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Figure  3.3:  Single  couples  used  to  represent  tensile  crack  with  dip  angle  8  and 
strike  angle  <pa  —  0,  (a)  (1,1)  with  M0  —  1,  (b)  (2,2)  with  M0  =  1,  (c)  (3,3)  with 
M0  = 

mation. 

The  potential  of  the  second  component  is  obtained  by  replacing  6  in  the 
potential  of  first  single  couple  with  6  +  90°.  The  third  potential  is  obtained  by 
rotating  the  first  single  couple  in  azimuthal  angle  90°  with  8  =  90°.  So  that,  the 
formulation  of  potential  for  tensile  crack  is 

3 

<f>  =  X^'(rA»  *;&»■&**) 

*=i 

Vv  =  ^2^'r{r,0k,z;8k,Mk)  (3.18) 

*=i 
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where  the  potentials  for  the  first  single  couple  (Fig  3.3(a)  are 
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and,  the  seismic  moment  Mo  =  A Au,  where  A  is  the  fault  area  and  u  is  the 
average  normal  displacement  to  the  fault.  The  formulation  for  a  tensile  crack 
with  an  oblique  strike  angle  is  obtained  by  rotating  the  azimuthal  angle  by  strike 
angle  <f>a  in  the  above  equation.  Since  the  particle  displacement  is  normal  to  the 
fault  surface  in  tensile  crack,  the  rake  angle  can  not  be  defined. 
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3.3  Global  Matrix  Approach 


Next  step  is  to  find  the  unknown  coefficients  in  Eq  3.6.  Denoting  the  set  of 
constraint  variables  for  each  Fourier  order  m  as 

wm(s,  z) 

um{s,z)  +  vm{s,z) 
um{s,z)  -vm{s,z) 

o?z{s,z)  +  o?z{s,z) 
o™{s,z)  -  ofc{s,z) 

the  boundary  condition  for  welded  contact  at  interface  n  separating  layers  n  and 
n  +  1  is 

IC (*,*)  +  K(s,zn)  -  K+1(s,0)  -P~  ,(8,0)  =  0,  (3.21) 

where  subscript  n  denotes  the  layer  number,  zn  is  the  thickness  of  layer  n,  and 
the  terms  with  tilde  ‘  ~  ’  denotes  source  terms.  Eq  3.21  gives  a  set  of  linear 
equations  for  each  interface.  These  sets  of  equations  for  each  interface,  when 
combined,  form  a  global  matrix  system,  which  is  solved  numerically.  To  be  more 
specific,  the  boundary  condition  can  be  rewritten  as 


An,iBn  A„+iiUBn+1  —  Rn+1)U  Rn,j5 


(3.22) 


where  the  coefficient  matrix  A  for  upper  interface  is 

—a  s  0  a  s  0 

—  5  (3  s  —s  —f3  s 

s  —(3  s  s  /3  s 

(2s2  —  k2)(i  —2 s/3^i  0  (2s2  —  k2)(j,  2s/3ft  0 

2 sap  -(2 s2  -  k2)n  —s(3fi  -2san  -(2s2  -  k2)n  s(3fi 

—2  sa/j,  (2s2  —  k2)n  —s/3fjL  2  sap  (2s2  —  k2)n  s/3fx 

(3.23) 
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and  the  matrix  for  the  lower  interface  is  obtained  by  multiplying  with  the  ap¬ 
propriate  exponential  matrix,  i.e. 


-A-n,j(s)  —  AniU(s)In(s) 


(3.24) 


where 


I n(s)  =  diag  {  e~a^Zn  ,  ,  e~^Zn  ,  ea^Zn  ,  e^z"  ,  e^Zn  }  .  (3.25) 


It  is  noted  that  the  coefficient  matrix  is  independent  of  Fourier  order  m,  which 
is  important  for  the  efficiency  of  the  numerical  code  since  the  multiple  right 
hand  side  due  to  Fourier  orders  in  global  matrix  can  be  solved  simultaneously, 
as  will  be  discussed  in  the  next  section  6.3.  The  unknown  variables  matrix  in 
the  potentials  is 


BrM  =  {<„M ,  CM  ,  CM  ,  CM  .  CM  .  CM  }  •  (3.26) 


The  source  matrices  R™  will  be  calculated  from  potentials  using  constitutive 
relations  for  displacements  and  stresses  for  different  types  of  sources.  The  source 
matrix  for  tensile  crack  is  given  below  and  other  source  matrices  are  included 
in  the  appendix.  The  following  source  matrices  are  found  from  the  potentials 
in  Eq  3.18  using  constitutive  relations  for  displacements  and  stresses,  which  are 
found  in  Eq  3.11  and  3.7.  For  cos  6  and  sin  26  Fourier  orders,  the  source  matrices 
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vanish.  For  order  m  =  0  , 


{0.25f(s2  -  2a2)  -  0.25f  cos  26* (s2  +  2a2)} 

+0.25f(l  +  3cos25fcs2e"^l1:-*'l 
0.25±  {(s2  -  2a2)  -  cos 2 6k(s2  +  2a2)} 

+0.25(1  +  3cos26k)f3se-^z~z^ 

—0.25^  {(s2  -  2a2)  -  cos  2Sk(s2  +  2a2)} 
r0  _  *  Mk  -0.25(1  +  3  cos  26k)0se-P\z~z-\ 

n  ~  hi  -0.25 {(s2  -  2a2)  -  cos  2 <5*(s2  +  a2)}  e~a\z~z *1 

— 0.5/is2/3(l  +  3cos24)e-^"2*l 
—0.5 (i$s  {(s2  —  2a2)  +  cos  2 Sk(s2  +  2a2}  e_“la'_z*l 
-0.25  /*fs(s2  +  /?2)(1  +  3  cos  26k)e-p\z~z-\ 

+0.5/ifs  {(s2  —  2a2)  —  cos  25*  (s2  +  2a2}  e~a I*-**! 
+0.25^s  (s2  +  /?2)(l  +  3  cos  28k)e-^z~z'\ 

(3.27) 


For  sinfl  order, 


ase-“lz-2*l  -  0.5  (s2  +  p2)je~p\z-z- 
fs2e~“  lz_z*l  —  fs2e~^z_z* 
Rl  _  ®  Mk  sin  2Sk  (  -^2e-“lz-z*l  +  02t~^z-z- 

jfc=2  4?rPnW2  —  /xfs(s2  +  /?2)e~“lz_z*l  +  /afs(s2  +  /?2)e_/?lz_z* 
— 2/^as2e-“lz-z*l  +  0.5//y  (3/?2  +  s2)e~^z_z* 
2/ias2e-“lz-z*l  -  0.5^(£  +  s2/?  +  2p3)e~^z-z- 


(3.28) 
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For  cos  29  order, 


p2  _  v'  Mt(l  -  cos 26k)  I 
n  “  4irpnu>2 


0.25cs2e-Q1lat-JB'l  -  0.25f*V 
0.25— e-alz-z*l  -  0.254e_^|z_2*1 

a  p 

— 0.25— e-alz_z*l  -  0.25^^se-^z-z*l 

a  P 

-0.25/i£(s2  +  /?2)e-“lz-**l  +  0.5/zs2/?e-^lz-2'l 
-0.5^sse~alz-z*l  +  0.5/xfs3e-/5lz-*'l 


0.5 /ifs3e  “lz  z*l  —  0.5 p$l32se 


'2se-e\*-*'  I 


In 


the  above  matrices,  M*  and  <5*.  are  defined  as 


(3.29) 


Mi  =  M0,  <5i  =  90°, 

M2  =  M0,  S2  =  6  +  90°,  and 

M3  =  ^M0,  4  =  ^, 

except  that  Mi  =  — M0  when  m  =  2,  cos  20  order.  The  by-product  of  the 
tensile  crack  formulation  is  an  explosive  source,  where  M*  =  Mo,  for  k  =  1,2, 3. 
Consequently,  only  the  zeroth  Fourier  order  component  is  non-zero,  yielding 
omnidirectional  radiation. 


The  local  matrix  in  Eq  3.21  for  each  interface  has  been  formulated  explicitly 
for  the  coefficient,  unknown  variables,  and  source  matrices  of  tensile  crack  as  an 
example.  Since  the  boundary  conditions  should  satisfied  at  all  the  interfaces, 
the  local  matrices  are  combined  to  form  a  global  matrix  (Fig  3.4).  Although  the 
dimension  of  the  global  matrix  is  increased  with  the  number  of  layers,  it  is  noted 
that  the  maximum  band  width  of  the  global  matrix  is  limited  to  constant  value 
18  so  that  a  special  band  matrix  solver  can  be  used  to  improve  the  numerical 
efficiency. 

Once  the  unknown  coefficients  are  found,  the  complete  solution  for  each  layer 
is  found  by  performing  inverse  Hankel  transform,  summing  over  Fourier  orders 
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Local  system 


Figure  3.4:  Mapping  between  local  and  global  systems  of  equations  (after 
Schmidt  and  Jensen,  1985). 


for  azimuthal  angle,  and  Fourier  synthesis  in  time  domain,  as  will  be  discussed 
in  the  next  section. 
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3.4  Field  Representation  in  a  Laterally  Strati¬ 
fied  medium 

3.4.1  Complete  Solution 

From  the  global  matrix  solution,  the  complete  solution  in  horizontal  wave  num¬ 
ber  domain  combining  homogeneous  and  particular  solution  is  found  for  each 
discretized  wave  number.  This  complete  solution  in  s-domain  is  called  Hankel 
transform  integrand  for  each  m-th  Fourier  order,  H™, 

H>,z,)  =  F"(«,2,)+P“(s,z,) 

=  An(s)Bn  (s,  zr)  +  It.,,  (s,  zr)  ,  (3.30) 

where  the  subscript  n  is  the  layer  number  in  which  the  receiver  at  z  =  zr  is 
located.  Note  that  the  integrand  is  not  a  function  of  r,  but  of  the  horizontal 
wave  number  s.  In  order  to  find  the  field  parameters,  inverse  Hankel  transform 
is  performed  on  H™(s,zr), 

H-(r,*r)  =  /  H-(s,*r)5Jm(r,s)d5 

Jo 


The  first  term  I  in  Eq  3.31  is  equivalent  to  the  terms  in  the  square  bracket  in 
Eq  3.7,  while  the  second  term  II  corresponds  to  Eq  3.27-  3.29  for  tensile  crack, 
except  that  they  are  now  evaluated  at  the  receiver  depth  z  —  zr  instead  of  at 
the  boundaries  z  =  zn.  Once  the  field  parameters  of  interest  for  Fourier  order 
m  in  Eq  3.31  are  found,  the  final  solution  is  found  from  Eq  3.5  summing  all  the 
Fourier  orders. 
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The  procedure  of  obtaining  the  complete  solution  in  real  space  1  has  been 
explained  in  detail.  Essentially,  the  solution  technique  is  quite  analytical,  until 
the  estimation  of  inverse  Hankel  transformation  is  introduced,  which  is  done 
numerically.  The  numerical  evaluation  of  the  integrals  is  important  to  the  ef¬ 
ficiency  of  the  developed  code.  Therefore,  the  aspects  of  numerical  integration 
will  be  discussed  in  the  next  section. 

3.4.2  Numerical  Consideration 

There  are  two  numerical  intensive  steps  in  the  solution  technique,  accompanied 
by  numerical  stability  and  efficiency  problems.  One  is  the  numerical  stability 
of  the  global  matrix  to  be  solved.  This  problem  is  treated  in  the  previous 
publication  [39],  and  the  numerical  stability  is  ensured  unconditionally  with 
proper  scaling  and  arrangement  of  the  coefficient  matrix  and  partial  pivoting 
when  solving  the  global  system  of  equation. 

The  other  is  the  inverse  Hankel  transform,  which  is  estimated  numerically. 
The  inverse  Hankel  transform  can  be  basically  estimated  in  two  ways,  FFT 
and  direct  numerical  integration.  The  FFT  method  is  an  approximation  to 
farfield  with  incident  angle  not  so  small,  since  the  inverse  Hankel  transform 
is  reduced  to  the  form  of  Fourier  transform  for  large  argument  of  the  Bessel 
function  [12]  [40].  The  advantage  of  this  method  is  the  computational  efficiency, 
since  the  reduced  Fourier  transform  is  computed  through  FFT  so  that  the  results 
are  found  at  every  range  step  simultaneously.  The  direct  integration  method, 
however,  requires  the  calculation  of  Bessel  function,  and  the  integration  scheme 
should  be  applied  for  each  argument  of  the  Bessel  function.  The  advantage  of 
integration  method  is  the  accurate  estimation  of  the  field  including  the  near 

1to  avoid  confusion  with  the  complete  solution  in  horizontal  wave  number  s-domain,  i.e. 
before  inverse  Hankel  transform. 
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field,  although  it  is  computationally  extensive.  Therefore,  this  method  can  be 
applied  efficiently  for  the  cases  of  a  few  receivers.  While  various  techniques  can 
be  used  to  increase  the  efficiency  of  the  numerical  integration  for  special  cases 
[2]  [11]  [19],  the  presently  available  methods  in  the  code  are  FFT  method  and 
direct  integration  methods  using  Simpson  rule  or  trapezoidal  rule  with  Romberg 
scheme. 


3.5  Numerical  Examples 

To  demonstrate  the  solution  technique,  the  radiation  in  a  homogeneous  un¬ 
bounded  medium  from  3  different  seismic  sources  are  considered.  These  exam¬ 
ples  are  chosen  due  to  the  simplicity,  and  the  existence  of  the  analytical  solution 
available  at  least  partially,  for  such  as  radiation  pattern. 

Case  3.1  Tensile  crack  :  The  tensile  crack  is  compact  and  located  at  the  origin 
with  receivers  at  2  =  60  m  (Fig  3.6),  giving  the  contour  level  of  vertical  displace¬ 
ment.  The  fault  orientation  parameters  are  6  =  90°  and  <f>t  =  0°.  The  source 
pulse  to  generate  synthetic  time  series  is  given  in  Fig  3.5.  Only  the  frequency 
components  between  0  -  300  Hz  are  considered. 

Case  3.2  Dip-slip  :  The  source  and  environmental  parameters  are  same  as  the 
case  3.2,  except  that  the  source  is  now  dip-slip. 

Case  3.3  Strike-slip  :  The  source  and  environmental  parameters  are  same  as 
the  case  3.2,  except  that  the  source  is  now  strike-slip. 

In  order  to  characterize  the  radiation  from  different  sources,  two  types  of 
displays  are  used.  One  is  the  contour  plot  of  a  field  parameter  at  certain  depth, 
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Modulus  (10  4)  Seismic  Moment  (N*m) 


SOURCE  PULSE 


Time  (seconds) 


SOURCE  SPECTRUM 


Frequency  (Hz) 


Figure  3.5:  Source  function  (a)  time  series,  (b)  spectral  shape. 


Compact  fault  surface  at  the  origin 


Figure  3.6:  Sketch  of  coordinate  system  with  source  and  receiver  positions. 

which  shows  the  radiation  pattern.  The  other  is  the  synthesized  time  series  for 
given  source  function.  These  are  discussed  subsequently. 

The  contour  plot  of  vertical  displacement  at  z  =  60  m  for  frequency  100  Hz 
corresponding  to  tensile  crack,  dip-slip,  and  strike-slip  sources  at  z  —  0  m  with 
dip  angle  6  =  90°  in  a  homogeneous  elastic  medium  are  given  in  Fig  3.7,  3.9, 
and  3.11  ,  respectively.  The  contour  level  is  taken  to  be 

CL  =  20log10  |u>|  re  1  Pa.  (3.32) 

The  source  strength  for  the  following  examples  is  the  unit  seismic  moment.  As 
expected  from  the  source  matrices,  the  Fourier  orders  excited  by  a  tensile  crack 
are  zeroth  and  second,  cos  29  orders  (Fig  3.7)  ,  while  only  first  order  is  excited 
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to  give  dipolar  radiation  pattern  in  horizontal  plane  for  dip-slip  (Fig  3.9)  and 
second  order  to  give  quadrupolar  radiation  pattern  for  strike-slip  (Fig  3.11)  . 
It  is  noted  that  the  radiation  patterns  for  different  seismic  sources  are  quite 
distinct  for  given  receiver  positions.  However,  the  strike-slip  can  be  obtained 
by  rotating  dip-slip  in  an  homogeneous  unbounded  medium,  which  suggests  the 
radiation  pattern  is  obviously  different  depending  on  the  receiver  positions. 


For  the  source  time  function  in  Fig  3.5,  the  synthetic  seismograms  for  dif¬ 
ferent  types  of  sources  are  given  in  Fig  3.8,  3.10,  and  3.12.  The  first  arrival  in 
Fig  3.8  for  tensile  crack  is  the  compressional  wave,  and  the  second  arrival  is  the 
shear  wave.  The  amplitude  of  the  compressional  wave  does  vary  slowly  while 
the  shear  wave  amplitude  vanishes  for  horizontal  angle  0  and  tt.  This  is  due  to 
the  fact  that  the  shear  wave  component  of  zeroth  and  second  orders  cancel  each 
other  at  these  directions.  The  seismogram  in  Fig  3.10  for  dip-slip  with  dip  angle 
8  =  90  consists  of  only  first  Fourier  order,  as  the  radiation  pattern  (Fig  3.9) 
suggests.  For  this  case,  the  shear  component  is  dominant  and  the  compressional 
component  is  barely  seen.  Also,  it  is  noted  that  the  phase  between  horizontal 
angle  0°  <  0  <  180°  are  shifted  by  180°  giving  negative  phase.  For  strike-slip 
with  dip  angle  8  =  90°  and  strike  angle  4>t  =  0°,  the  second  Fourier  order  is 
the  only  contribution  as  in  Fig  3.11.  Again  phases  are  alternating  every  90°. 
Therefore,  it  is  argued  that  the  temporal  characteristics,  as  well  as  radiation 
patterns,  of  different  sources  can  be  used  to  identify  the  source  parameters. 

For  above  three  cases  in  a  homogeneous  unbounded  medium,  the  compres- 


(3.33) 

(3.34) 
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Range  (km) 


Figure  3.7:  Radiation  pattern  of  tensile  crack  with  dip  angle  S  =  90°  for  single 
frequency  /  =  100  Hz  in  a  homogeneous  elastic  medium. 


Time  (seconds) 

Figure  3.8:  Seismogram  for  tensile  crack  with  dip  angle  6  =  90°,  receivers  at 
z  =  0  m,r  =  300  m  at  every  30°  in  a  homogeneous  elastic  medium. 
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Figure  3.9:  Radiation  pattern  of  dip-slip  with  dip  angle  S  =  90°  for  single 
frequency  /  =  100  Hz  in  a  homogeneous  elastic  medium. 
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Figure  3.10:  Seismogram  for  dip-slip  with  dip  angle  6  =  90°,  receivers  at 
z  —  0  m,  r  =  300  m  at  every  30°  in  a  homogeneous  elastic  medium. 
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Range  (km) 


Figure  3.11:  Radiation  pattern  of  strike-slip  with  dip  angle  6  =  90°  for  single 
frequency  /  =  100  Hz  in  a  homogeneous  elastic  medium. 


Time  (seconds) 

Figure  3.12:  Seismogram  for  strike-slip  with  dip  angle  6  =  90°,  receivers  at 
z  —  0  m,r  —  300  m  at  every  30°  in  a  homogeneous  elastic  medium. 
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where  rr  and  zr  are  the  coordinates  of  receiver  location.  For  the  cases  considered, 
since  the  medium  is  unbounded,  there  is  no  dispersion  of  the  existing  waves. 
However,  when  there  wave  guides  are  introduced  by  boundaries  as  in  the  case 
of  a  floating  ice  plate  in  a  central  Arctic  environment,  the  group  velocity  is 
different  from  the  phase  speed  of  corresponding  waves.  These  cases  are  discussed 
in  Chapter  6. 
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Chapter  4 

Radiation  from  Propagating 
Cracks 


The  excitation  by  a  seismic  source  originates  on  the  dislocation  surface.  Often, 
the  size  of  dislocation  surface  is  approximately  larger  than  a  quarter  of  the  wave 
length  of  interested  frequency,  causing  directivity  pattern,  as  has  been  observed 
in  some  earthquakes  [17]  .  The  directivity  pattern  is  caused  by  the  phase  in¬ 
terference  introduced  by  the  geometry  of  crack  surface,  and  crack  propagation 
speed  or,  rupture  speed  in  seismology.  Obviously,  when  the  seismic  fault  dimen¬ 
sion  is  relatively  small  enough  not  to  cause  the  phase  interference,  they  can  be 
treated  as  point  or  compact  sources. 

In  Chapter  4.1,  the  solution  to  a  canonical  problem  of  non-compact  crack 
radiation  shows  that  the  formulation  of  radiation  pattern  is  actually  the  same 
as  that  of  an  array  of  sonic  transducers.  Therefore,  the  developed  model  can 
be  also  applied  to  the  radiation  from  an  array  of  sources  as  well  as  propagating 
cracks.  In  Chapter  4.2,  an  efficient  numerical  model  is  developed  for  the  radia¬ 
tion  and  propagation  from  propagating  cracks  in  a  laterally  stratified  medium. 
In  the  numerical  examples  section,  the  same  canonical  problem  is  solved  using 
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L :  Gack  length 


►  x 


Figure  4.1:  Discretized  source  distributed  over  the  line  crack  and  coordinate 
system. 

the  developed  model,  and  the  characteristics  of  radiated  field  from  propagating 
cracks  are  discussed. 


4.1  Analytic  Solution  to  a  Canonical  Crack  Ra¬ 
diation  Problem 

A  canonical  problem  simulating  the  radiated  field  in  an  homogeneous  unbounded 
acoustic  medium  can  be  analytically  solved,  giving  insights  for  the  radiation  in 
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a  laterally  stratified  elastic  medium.  The  fault  surface  1  can  be  modeled  as  an 
array  of  discretized  sound  sources  as  in  Fig  4.1,  where  a  line  crack  is  considered. 
The  total  pressure  field  with  contributions  from  the  sources  can  be  expressed  as 


P(x>  t) 


AAl  r. 

it  r(i 


,  Al  A  l'  r.  .  A  l  A  V 
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+  •••  1  , 
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c 
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(4.1) 


where  w,-(t)  is  the  source  function  of  the  t-th  discrete  source,  L  is  the  crack  length 
(array  length) ,  A  is  the  source  strength  of  an  array,  cv  is  the  crack  propagation 
speed  (phase  velocity  for  an  array),  Al  is  the  distance  between  the  discretized 
sources,  and  Al'  is  the  distance  between  adjacent  discrete  sources  to  the  receiver 
(A l'  =  A/sinxcos(0  —  $<*)).  Assuming  that  each  discrete  source  has  the  same 
source  function,  the  pressure  field  can  be  written  as 


,  .  AAl  ”  ..  nAL  ) 

p(x>*)  =  7—  £  «>((*-—- )--  „ 

^ _ KT  C 


nAl' 


)• 


(4.2) 


n——N 


As  can  be  seen  in  the  above  equation  4.2,  there  are  two  phase  terms  introduced. 
The  first  term  nAl/cv  is  the  time  delay  taken  for  the  crack  tip  to  travel  the 
distance  between  discretized  sources  nAl.  The  second  term  —nAl'/c  is  the 
time  delay  due  to  the  different  source  location  to  the  receiver  positions.  For 
the  limiting  case  of  Al  — *■  0,  i.e.  for  the  continuously  distributed  sources,  the 
expression  in  Eq  4.2  reduces  to  integral  form 


,  ,  A 1  [L /2  ,  .  1 

P(x,0  =  7- /  w{t-ra{- 

JLj  T  Cfj 


sinxco .{>-)<)  _  r  ^ 
c  c 


1Even  if  there  can  not  exist  a  fault  surface  in  an  acoustic  medium,  the  fault  surface  denotes 
the  general  surface  over  which  the  sources  are  distributed. 
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Defining 


(4.4) 


„  1  sin  x  cos  (0  —  6d) 

cv  c  ’ 

the  pressure  expression  can  be  written  as 

A  1  [L/2  f 

P(X’()  =  Jr  J-L/l  3f“  ~~  <4  ir‘ 

=  4-  /  <"(<  -  Hr.  -  -)  rect(4)*a  ,  (4.5) 

Li  r  J  OO  C  Z 

where  the  rectangular  function  rect(y)  is  defined  as 

1  ,  when  |r0|  <  f 
0  ,  otherwise 

Using  transformation,  E ra  —  r, 

,  .  ^41  fLl 2  .  r.  1  ,L.  , 

p(x,t)  =  v-/  u>(t  —  r - )  — rectf— )  dr 

'  LrJ-L/2  v  c  E  v2' 

A1  .  r.  1  .L. 

=  *  §rect(2^  ^ 
The  rectangular  box  is  the  result  of  uniform  source  strength  distribution  along 
the  crack  surface  using  farfield  approximation.  Of  course,  the  source  strength 
function  does  not  have  to  be  a  box  type,  it  can  be  a  triangle,  or  any  kind  of  shape 
that  best  describe  the  source  strength  distribution.  Since  the  source  distribution 
function  is  not  known,  the  rectangular  shape  is  assumed  for  simplicity. 

From  Eq  4.6,  it  is  seen  that  the  signal  observed  at  the  hydrophone  is  the 
source  function  convolved  with  the  source  distribution  function,  which  is  a  func¬ 
tion  of  E  =  E(0d,x,c«j)  as  shown  in  Eq  4.4,  so  that  the  observed  signal  will  be 
different  for  different  receiver  directions.  However,  the  area  under  the  rectangu¬ 
lar  function  remains  constant  for  all  the  receiver  directions. 

Since  the  convolution  in  time  domain  corresponds  to  the  multiplication  of 
two  Fourier  transforms  in  the  frequency  domain,  it  will  be  useful  to  look  at 
the  frequency  domain  solution.  In  order  to  find  the  frequency  domain  solution, 


rect(-)  = 
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the  impulse  source  function  is  assumed  i.e.  w(t)  =  6(t).  The  plot  of  frequency 
domain  solution  is  given  in  Fig  4.2  (b).  The  Fourier  transform  of  a  rectangular 
function  is  sine  function  (sinc(r)  =  sinr/r),  so  that  the  peaks  and  nulls  as 
function  of  r  corresponds  to  the  radiation  pattern  in  the  real  space  for  certain 
range  of  r,  which  is  determined  by  the  crack  propagation  speed  and  the  ratio  of 
wave  length  to  the  crack  dimension.  In  fact,  the  formulation  is  exactly  same  as 
the  radiation  from  the  sonar,  which  makes  the  model  applicable  to  the  radiation 
from  an  array  of  transducers,  such  as  airguns  or  explosions. 

The  interpretation  of  the  results  for  this  canonical  problem  gives  useful  in¬ 
formation  about  the  spectral  and  temporal  radiation  characteristics.  First,  the 
frequency  domain  solution  shows  that  the  spectral  level  at  fixed  receiver  posi¬ 
tion  (  i.e.  fixed  5)  is  different  for  different  frequencies,  which  is  called  frequency 
dependent  directivity.  This  is  one  of  the  characteristics  that  distinguishes  the 
non-compact  sources  from  compact  sources.  The  time  domain  solution  in  Eq  4.6 
shows  the  convolution  of  source  signal  with  different  duration  of  rectangular  box 
depending  on  the  receiver  directions.  These  aspects  of  the  frequency  and  time 
domain  solutions  are  discussed  in  the  numerical  examples  Case  4.1,  where  the 
same  example  is  treated  for  a  given  source  function  using  the  presently  developed 
code. 


4.2  Numerical  Model  for  Propagating  Cracks 
in  an  Stratified  Medium 

The  field  caused  by  the  moving  crack  can  be  basically  found  by  integrating 
the  proper  Green’s  function  of  the  specific  sources  over  the  fault  surface.  This 
method  assumes  that  there  is  no  scattering  from  the  pre-existing  crack  surface, 
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Figure  4.2:  Fourier  transform  pair,  (a)  rectangular  source  strength  distribution 
function  along  the  crack,  (b)  Fourier  transform  of  ^rect(ra). 
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Figure  4.3:  Sketch  of  moving  crack. 

i.e.  the  distributed  sources  over  the  crack  surface  radiate  as  if  the  medium  is 
homogeneous. 

In  order  to  treat  the  moving  dislocation  sources,  it  is  necessary  to  place  the 
sources  off  the  2- axis  in  Fig  4.3  ,  which  has  been  studied  in  earlier  development 
of  the  three  dimensional  version  of  SAFARI  code  [38]  .  However,  the  mathemat¬ 
ical  treatment  is  rather  complex  and  requires  the  calculation  of  higher  orders  of 
Bessel  functions  causing  convergence  problem  and  numerical  inefficiency.  There¬ 
fore,  simple  transformation  of  source  and  receiver  position  is  used  in  this  study. 
The  actual  problem  of  distributed  sources  over  the  fault  surface  is  reduced  to 
a  numerical  problem  with  one  source  on  2-axis,  and  corresponding  receivers  for 
each  of  the  sources  (Fig  4.4). 

The  total  radiation  field  due  to  moving  crack  can  be  obtained  by  summing  the 
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b)  reduced  numerical  model. 


Figure  4.4:  Receiver-source  position  transformation  (a)  source  distribution  over 
fault  surface,  (b)  reduced  numerical  model. 
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fields  at  receiver  positions.  In  order  to  apply  transform-superposition  method, 
some  considerations  are  needed.  First,  the  Green’s  function  of  the  source  po¬ 
sitioned  on  the  z-  axis  has  to  be  multiplied  by  a  phase  term  due  to  traveling 
crack  speed.  The  phase  difference  introduced  by  a  propagating  crack  consists  of 
two  parts.  One  corresponds  to  the  time  delay  for  the  crack  to  propagate  from 
one  source  position  to  another,  and  the  other  is  a  time  delay  due  to  different 
source  positions,  as  discussed  in  the  previous  section.  Therefore,  the  total  phase 
difference,  for  example  (Fig  4.3),  caused  by  different  source  positions,  rj  and  r2, 
is 


To 


L  ,  kzl-kil 

vr  c 


(4.7) 


The  second  term  is  included  in  receiver-source  position  transformation.  The  first 
term  is  a  phase  to  be  multiplied  to  the  source  Green’s  function.  Secondly,  when 
the  field  parameters  are  added,  another  transformation  to  correct  the  receiver- 
source  transformation  is  necessary.  This  consideration  is  made  to  resolve  the 
farfield  approximation  for  the  array,  i.e.  the  following  requirement  for  farfield 
approximation  is  not  necessary 


L  <C  Vr2  +  z2  , 


(4.8) 


of  which  restriction  applies  to  the  analytical  solutions. 

For  the  case  of  the  multiple  depth  sources  (Fig  4.5),  the  same  procedure 
will  be  performed  for  each  source  depth.  This  is  due  to  the  fact  that  the  ge¬ 
ometry  concerned  is  a  laterally  stratified  medium,  so  that  the  source-receiver 
transformation  can  be  applied  to  horizontal  direction  only. 
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Figure  4.5:  Two  dimensional  crack  surface  and  source  distribution. 

4.3  Numerical  Examples 

Two  canonical  problems  are  considered  in  acoustic  and  elastic  media.  The  media 
are  homogeneous  and  unbounded,  so  that  there  is  no  boundary  effect.  This  will 
allow  us  to  look  at  the  effects  of  non-compact  sources.  The  source  and  receivers 
of  examples  are  explained  in  the  followings. 

Case  4.1  Radiation  from  an  array  in  a  homogeneous  unbounded  acoustic  medium 
:  The  line  array  is  extended  from  (—25,0,0)  to  (25,0,0)  m,  traveling  at  speed 
cv  =  1800  m/sec.,  while  the  sound  speed  in  the  acoustic  medium  is  c  = 
1440  m/sec.  The  source  type  is  explosive  source  (monopole).  Refer  to  Fig  4.3. 
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Case  4.2  Radiation  from  a  propagating  crack  in  a  homogeneous  unbounded  elas¬ 
tic  medium  :  The  crack  is  assumed  to  be  a  line  crack  with  a  element  Green’s 
function  being  a  explosive  source,  extending  from  (—25,0,0)  to  (25,0,0)  m.  The 
crack  travels  at  cv  =  1800  m/sec,  and  the  compressional  and  shear  speed  of  the 
medium  are  ce  —  3500  m/sec,  c,  =  1800  m/sec,  respectively. 

Two  types  of  displays  are  used  to  visualize  the  radiated  field.  For  Case  4.1, 
the  contour  plot  for  xy— plane  at  z  =  60  m  (  Refer  to  Fig  3.6  )  and  the  synthetic 
time  series  at  range  300  m  are  given  in  Fig  4.6  and  4.7,  respectively.  In  Fig  4.6, 
the  major  lobe  direction  0m  due  to  steering,  that  is  due  to  propagation  speed, 
can  be  calculated  simply  by  Snell’s  law, 


This  major  lobe  corresponds  to  the  peak  of  Sine  function  in  Fig  4.2  shifted  by  a 
corresponding  phase  due  to  crack  propagation  speed,  or  phase  speed  of  an  array. 
The  synthetic  time  series  in  Fig  4.7  shows  the  major  lobe  effect  in  the  direction 
at  «  ±60°  on  arrivals  with  high  amplitude.  When  the  higher  frequencies  with 
the  same  crack  dimension  are  considered,  the  steering  effect  will  become  more 
obvious  radiating  the  acoustic  energy  spatially  in  a  narrower  beam. 

For  Case  4.2,  the  elastic  medium  has  two  wave  components,  compressional 
and  shear  waves.  However,  since  the  source  is  omni-directional  explosive  source, 
the  shear  wave  is  not  excited,  which  simplies  the  problem  2.  The  compressional 
wave  speed  is  3500  m/sec,  which  is  greater  than  the  crack  propagation  speed 
1800  m/sec,  called  subsonic  crack  speed.  For  this  case,  the  steering  angle,  or 
the  major  lobe  direction  0m  can  not  defined  in  real  space,  i.e.  0m  is  imaginary 

2  The  shear  wave  excited  by  directional  sources  is  treated  in  Chapter  6 
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Range  (km) 


Figure  4.8:  Radiation  pattern  for  an  array  of  explosive  sources  for  a  single 
frequency  f  =  100  Hz  in  a  homogeneous  elastic  medium. 


Time  (seconds) 

Figure  4.9:  Seismogram  of  the  vertical  velocity  for  an  array  of  explosive  sources 
with  receivers  at  z  —  60  m,  r  =  300  m  at  every  30°  in  a  homogeneous  elastic 
medium. 
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in  Eq  4.9.  Therefore,  the  major  lobe  is  not  seen  in  Fig  4.8,  and  the  directivity 
pattern  is  not  as  drastic  as  Case  4.1  since  the  ratio  of  the  compressional  wave 
length  to  crack  length  is  greater  than  the  previous  case.  This  effect  is  shown 
in  Fig  4.9  with  moderate  changes  in  amplitude  between  channels  compared  to 
Fig  4.7,  which  means  the  beam  is  broader  spatially. 
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Chapter  5 

Transversely  Isotropic  Medium 


In  this  Chapter,  it  is  shown  how  the  global  matrix  method  can  be  applied  to  the 
formulation  of  propagation  in  a  transversely  isotropic  medium.  A  transversely 
isotropic  medium  is  characterized  by  five  independent  elastic  constants.  Often, 
the  situation  occurs  when  the  periodically  repeated  finely  layered  sediment,  or 
the  sea  and  lake  ice  are  considered. 

Since  the  theory  for  the  wave  propagation  in  a  transversely  isotropic  medium 
is  well  established  [3]  [13]  [14]  [15]  [20]  [26]  [27]  [28],  the  emphasis  is  given  on  the 
application  of  global  matrix  method  to  a  transversely  isotropic  medium.  First, 
the  equation  of  motion  is  presented  in  cylindrical  coordinates.  The  equation  of 
motion  is  no  longer  expressed  as  a  wave  equation,  as  in  the  case  of  an  isotropic 
medium.  Then,  the  equivalent  transversely  isotropic  elastic  constants  for  pe¬ 
riodically  fine  layered  medium  are  summarized.  By  introducing  intermediate 
functions,  the  decoupled  equations  are  formulated  for  SH  and  SV  —  P  waves. 
These  functions  are  utilized  to  express  the  boundary  conditions  at  the  interfaces 
in  a  form  compatible  with  the  isotropic  case.  It  is  noted  that,  in  order  to  set  up 
a  system  of  linear  equations,  the  intermediate  functions  need  to  be  arranged  as 
the  potentials  in  an  isotropic  medium,  which  is  discussed  in  Section  5.4. 
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5.1  Equation  of  Motion 


The  equation  of  motion  in  cylindrical  coordinates  is  [43,  pp.219] 

P-QP  =  pF  +  (div  oT  ,  div  o9  ,  div  oz  )  +  -(- o H  ,  or+  ,  0  )  .  (5.1) 

For  an  isotropic  medium,  the  above  equation  5.1  reduces  to  Eq  2.8  using  the 
stress-strain  relationship  with  two  independent  elastic  constants  A  and  n.  In 
the  following  sections,  the  general  equation  of  motion  5.1  will  be  further  sim- 
plied  for  a  transversely  isotropic  medium  to  a  system  of  equations  with  the  five 
independent  elastic  constants,  which  decouples  for  P  —  SV  and  SH  waves  by 
introducing  a  set  of  proper  intermediate  functions,  f. 


5.2  Stress-Strain  Relation 

The  stress-strain  relation  in  cylindrical  coordinates  is 
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where  C12  =  Cu  —  2C66.  The  expressions  for  strain  in  cylindrical  coordinates 
are 
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Figure  5.1:  Sketch  of  a  periodically  finely  layered  medium. 

of  which  general  relation  in  curvilinear  coordinates  are  found  in  Takeuchi  and 
Saito  [43].  For  an  isotropic  medium,  the  stress-strain  relation  can  be  found  by 
substituting  Cn  =  C33  =  A  +  2 fj,,  Cu  =  C66  =  n,  and  C12  =  C1S  =  A. 

5.3  Equivalent  Transversely  Isotropic  Elastic  Con¬ 
stants  for  Periodically  Fine  Layered  Medium 

The  periodic  finely  layered  medium  described  in  Fig  5.1,  when  H  is  much  smaller 
than  the  wave  length  of  the  interested  frequency,  can  be  treated  as  an  equiva¬ 
lent  homogeneous  transversely  isotropic  medium.  Consequently,  the  equivalent 
elastic  constants  for  a  periodically  finely  layered  medium  with  N  homogeneous 
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layers  (Fig  5.1)  may  be  shown  to  be 


C44  =  (M-1)-1 

C66  =  (m) 

C33  =  {' l/l*)~ 1 

C13  =  (l-2{7))(7/M)-1 

<?ii  =  4(/i)  -4('tn)  +  (l-2(7))2(7/^)-1 

C12  =  Cxi  —  2{*i)  =  Cu  —  2C66  , 

where  the  angled  bracket  denotes  the  thickness- weighted  average  [41]. 


(5.4) 


5.4  Decoupled  Equations  for  SH  and  P-SV  Waves 

Returning  to  solving  the  equation  of  motion,  Eq  5.1  can  be  reduced  to  the  equa¬ 
tion  with  respect  to  displacement  u,  first  using  stress-strain  relation  and  then, 
using  strain-displacement  relation  as  in  Eq  5.2.  While  the  three  displacement 
potentials,  <j>.  A,  and  rp,  defined  in  Eq  3.6  and  3.4,  that  express  the  P,  SV,  and 
SH  waves  respectively,  can  be  defined  in  an  isotropic  medium,  the  potentials 
cannot  be  defined  for  the  equation  of  motion  in  an  transversely  isotropic  medium 
since  the  P  —  SV  waves  cannot  be  separated.  However,  a  set  of  intermediate 
functions,  fi,  ft,  f 3,  and  /4,  which  combines  the  P  —  SV  waves  and  /6  and  /6  for 
the  SH  wave,  may  be  defined,  since  the  P  —  SV  and  SH  waves  still  decouple. 
Therefore,  the  problems  for  P  —  SV  and  SH  waves  are  solved  separately. 

In  cylindrical  coordinates,  the  solution  to  the  propagation  of  P  —  SV  waves 
may  be  assumed  as 

w  =  fz(z;u>,s)Y™  eiut 
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(5.5) 
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For  the  SH  wave,  the  solution  assumes  the  following  forms, 
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where  Y™  in  f$  and  /6  is  defined  as 

sin  mO 
—  cos  mO 

Comparing  Eq  5.8  with  5.6,  the  different  arrangements  of  Fourier  orders  for  the 
SV  —  P  and  SH  waves  are  similar  to  those  defined  in  Eq  3.4  in  an  isotropic 
medium  x.  Obviously,  this  is  to  enable  the  boundary  conditions  to  be  compat¬ 
ible  with  the  previous  formulation  in  the  isotropic  medium.  Substituting  these 

1The  major  difference  is  that  the  potentials  for  the  isotropic  medium  represent  the  P,  SV, 
and  SH  waves  fully  separated,  while  the  intermediate  functions  for  the  transversely  isotropic 
medium  represent  the  coupled  P  —  SV  waves  and  decoupled  SH  wave. 


(5.8) 
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equations  into  the  equations  of  motion,  it  can  be  shown  that 
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(5.9) 


The  above  equations  are  taken  from  the  references  [43]  [29],  and  rearranged  to 
be  compatible  with  global  matrix  method.  As  can  be  seen  in  the  above  Eq  5.9, 
/,-,  for  i  =  1,  ...4  are  decoupled  from  /,-,  for  i  =  5,6.  The  eigenvalues  for  the  first 
4  equations  can  be  found  from 
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(5.10) 

The  four  eigenvalues  represent  the  up-going  and  down-going  components  of 
quasi-P  and  quasi-SV'  waves.  The  two  eigenvalues  obtained  from 


i/2  = 


Nk2 


pu * 


(5.11) 


represent  the  up-going  and  down-going  components  of  the  SH  wave.  In  the 
isotropic  case,  ux  and  v2  corresponds  to  a  and  /?,  which  are  the  vertical  wave 
numbers  for  compressional  and  shear  waves,  respectively.  And,  U3  corresponds 
to  /3  since  the  SH  wave  has  the  same  wave  number  as  the  S'V'  wave  in  an 
isotropic  medium. 

In  Eq  5.10  and  5.11,  considering  the  vertical  wave  number  v  —*■  0,  in  which 
case  the  waves  travel  in  the  horizontal  direction,  the  phase  speeds  are  found  to 
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be  JCuJp  for  P  wave  and  \fC44Jp  for  SV  wave,  and  Jc^/p  for  SH  wave, 
respectively.  Now,  considering  k  — *  0,  in  which  the  waves  propagate  in  the 
vertical  direction,  the  phase  speeds  are  C33 / p  for  the  compressional  wave  P, 
and  yJCu/p  for  the  shear  waves,  SV  and  SH. 

Once,  the  eigenvalues  and  eigenvectors  are  found,  the  homogeneous  solution 
to  Eq  5.9  can  be  written 

f  =  EAK  ,  (5.12) 

where  E  is  the  eigenvector  matrix  of  Eq  5.9,  A  is  a  diagonal  matrix  for  eigenval¬ 
ues,  and  K  is  a  unknown  coefficient  matrix  for  six  components  of  up-going  and 
down-going  P,  SV,  and  SH  waves,  respectively.  They  are  written  as  follows, 
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where  7,-,  X,-,  and  Y,  are  defined  as 

_  ki/j(Ci3  +  C44) 

pu2  -  k2Cu  +  i/fCzz 

Xi  =  CzzUi'ii  —  kCiz  (5.16) 

Yi  =  Cu(ui  +  k'li)  . 

It  is  noted  that  A,  K  and  E  correspond  to  I,  B  and  A  for  the  isotropic  case  in 
Eq  3.25,  3.26,  and  3.23,  respectively.  As  in  the  case  in  the  isotropic  medium, 
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the  E  matrix  is  independent  of  Fourier  order  m,  while  the  unknown  coefficient 
matrix  K  has  to  be  found  for  each  orders  of  m.  However,  note  that  the  f  matrix 
is  an  intermediate  set  of  functions  like  potentials  in  the  isotropic  medium.  In  the 
next  Section,  the  intermediate  functions  are  combined  to  represent  the  boundary 
conditions  used  in  the  interface  in  the  global  matrix  method. 


5.5  Formulation  of  Boundary  Conditions  Com¬ 
patible  with  Global  Matrix  Method 


The  boundary  conditions  at  the  interfaces  (Eq  3.20)  are  formulated  based  on 
Eq  5.12,  which  is  a  form  of  homogeneous  solution  in  a  transversely  isotropic 
medium. 
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where  the  explicit  expression  for  f  can  be  calculated  from  Eq  5.12 
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Figure  5.2:  A  single  force  as  a  source. 

In  the  implementation  of  the  numerical  code,  the  eigenvectors  are  scaled  in  the 
same  way  as  the  coefficient  matrix  A  in  Eq  3.23  for  the  isotropic  medium  for 
compatibility. 


5.6  Numerical  Examples 

One  of  the  most  characteristic  features  of  the  wave  propagation  in  transversely 
isotropic  medium  is  the  separation  of  the  SH  and  SV  waves  depending  on 
the  direction  of  propagation.  The  separation  is  maximized  generally  for  the 
propagation  in  the  horizontal  direction,  as  can  be  seen  in  the  slowness  diagram 
in  Fig  5.3  which  has  been  computed  for  the  Case  5.1.  In  order  to  see  the 
separation,  the  source  needs  to  excite  the  SH  wave,  of  which  sources  include 
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Figure  5.3:  Slowness  surface  of  P(solid  line),  SV (dotted  line),  and  SH (dashed 
line)  waves  for  case  1. 

forces,  or  couples  with  Fourier  orders  greater  than  the  zero-th.  For  simplicity, 
a  single  force  in  an  unbounded  and  homogeneous  medium  2  is  considered  in  the 
following  examples. 

Case  5.1  Transversely  isotropic  medium  :  The  transversely  isotropic  medium 
considered  consists  of  two  isotropic  layers,  which  are  repeatedly  fine  layered. 
For  the  first  and  second  layer,  the  compressional  wave  speeds  are  4000  and 
2000  m/sec  ,  and  the  shear  speeds  2000  and  1000  m/sec  ,  respectively.  The 
density  are  1  g/cms  for  both  layers.  The  point  force  is  placed  at  the  origin  and 
the  direction  of  the  force  is  to  the  horizontal  angle  30°,  and  the  vertical  angle  is 

2  Since  the  presently  developed  codes  can  not  have  a  source  in  the  anisotropic  medium,  an 
infinitesimally  thin  isotropic  layer  is  placed  at  the  depth  of  source. 
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60°  towards  the  2-axis  as  shown  in  Fig  5.2.  The  12  receivers  are  placed  at  every 
30°  at  depth  of  60  m  and  range  of  300  m,  as  in  Fig  3.6. 

Case  5.2  Isotropic  medium  :  In  this  example,  the  isotropic  medium  is  consid¬ 
ered.  The  other  parameters  are  the  same  as  previous  Case  5.1,  except  that  the 
compressional  and  shear  wave  speeds  are  3000  and  1500  m/ sec,  respectively. 

Through  these  examples,  the  propagation  effects  in  a  transversely  isotropic 
medium  are  discussed.  The  slowness  surfaces  show  the  varying  phase  speeds 
depending  on  the  propagation  direction.  For  the  horizontal  direction,  the  waves, 
P,  SV,  and  SH,  travel  with  different  speeds,  while  for  the  vertical  direction, 
the  shear  waves,  SV  and  SH  waves,  travels  at  the  same  speed  so  that  the  shear 
waves  does  not  separate.  This  can  be  demonstrated  in  the  synthetic  time  series. 
Again,  the  same  source  function  in  Fig  3.5  is  used.  The  synthetic  time  series 
of  the  velocities  in  z,x,  and  y-directions  are  shown  in  Fig  5.4  and  5.5  for  the 
cases  5.1  and  5.2,  respectively.  The  comparison  of  the  seismogram  for  Case  5.1 
in  Fig  5.4  with  Case  5.2  in  Fig  5.5  shows  the  separation  of  the  SH  and  SV 
waves,  as  expected.  The  separation  time  can  be  calculated  using  the  receiver 
directions  from  the  slowness  surfaces  for  SV  and  SH  waves  in  Fig  5.3.  Since 
the  receiver  positions  are  almost  horizontal,  SH  arrives  faster  than  SV  wave. 
The  relative  strength  of  each  modes  of  propagation  as  well  as  the  propagation 
speed  in  the  horizontal  direction  can  be  also  identified  in  the  modulus  plot  of 
the  Hankel  inverse  transform  integrand  in  Fig  5.6. 

First,  the  equivalent  five  elastic  constants  can  be  calculated  for  the  case  5.1 
from  Eq  5.4. 

C44  =  1.6  x  109  ,  C66  =  2.5  x  109  , 

C33  =  6.4  x  109  ,  C1S  =  3.2  x  109  ,  (5.19) 

Cn  =  9.1  x  109  , 
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Figure  5.5:  Synthetic  time  series  in  an  unbounded  homogeneous  isotropic 
medium,  (a)  vertical  velocity,  (b)  x-direction  velocity. 
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Figure  5.6:  Modulus  of  the  inverse  Hankel  transform  integrand  for  the  isotropic 
medium  (solid  line),  and  transversely  isotropic  medium  (broken  line),  (a)  and 
(b)  show  the  vertical  displacements  for  Fourier  orders  zero-th  and  cos  6  orders, 
respectively,  (c)  and  (d)  for  the  x-direction  velocity,  and  (e)  and  (f)  for  the 
y-direction  velocity. 
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where  the  unit  is  Kg /m- sec2.  From  these  elastic  constants  and  using  the  center 
frequency  of  100  Hz,  the  phase  speeds  and  wave  numbers  for  the  P,  SV,  and 
SH  waves  in  the  horizontal  direction  are  found  as 
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For  case  5.2,  the  phase  speeds  for  all  wave  types  are  the  same  in  all  directions, 
and  the  horizontal  wave  number  can  be  calculated  from  the  phase  speed  and 
frequency, 
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(5.21) 


Now,  the  moduli  of  the  inverse  Hankel  transform  integrand  in  Fig  5.6  can  be 
interpreted  based  on  the  horizontal  wave  numbers  found  from  the  phase  speed  of 
each  mode.  The  peaks  at  s  =  0.419  and  0.474  in  Fig  5.6(a)  for  the  zeroth  Fourier 
order  of  vertical  velocity  represent  the  dominant  SV  wave  in  the  isotropic  and 
the  transversely  isotropic  medium,  respectively.  For  the  first  Fourier  order  of 
the  vertical  velocity  in  Fig  5.6  (b),  the  peak  at  s  =  0.474  for  the  SV  wave  in  the 
transversely  isotropic  medium  is  observed,  which  does  not  exist  in  the  isotropic 
medium.  The  magnitude  of  the  zeroth  order  is  greater  than  the  first  order  by 
an  order  of  10,  the  directional  variation  in  the  horizontal  angle  is  negligible  in 
the  synthetic  time  series  for  the  vertical  velocity  in  Fig  5.4  and  5.5  (b)  and  (c). 
The  peaks  at  s  =  0.200  and  0.209  in  Fig  5.6  (d)  for  the  first  Fourier  order  shows 
the  P  wave  in  the  isotropic  and  the  transversely  isotropic  medium,  respectively. 
Another  peak  at  5  =  0.379  is  the  contribution  of  the  SH  wave,  of  which  peak 
appears  only  in  the  transversely  isotropic  medium. 
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Although  the  modulus  of  the  inverse  Hankel  transform  integrand  is  not  the 
final  form  of  the  quantitative  field  presentation,  it  gives  the  idea  of  what  modes 
of  wave  exist  and  are  dominant.  The  previously  given  examples  are,  however, 
the  most  simple  cases  without  any  boundaries,  where  the  group  velocity  is  the 
same  as  the  phase  speed,  and  not  dispersive.  When  there  exists  a  wave  guide,  the 
Hankel  transform  integrand  can  be  very  useful  in  identifying  the  phase  speed  of 
existing  modes  along  with  the  dispersion  relation.  These  examples  in  relation  to 
the  wave  propagation  in  the  floating  ice  plate  in  the  central  Arctic  environment 
will  be  discussed  in  Chapter  6.5. 
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Chapter  6 

Application  to  Ice  Crack 
Radiation  in  the  Central  Arctic 
Environment 


The  main  source  of  the  central  Arctic  ambient  noise  has  been  known  to  originate 
from  the  elastic  motion  of  the  ice  plate  caused  by  environmental  stresses.  When 
the  stress  in  the  ice  plate  is  locally  greater  than  the  strength  of  the  ice  causing 
fracture,  the  stored  energy  is  released  in  the  form  of  elastic  motion  of  the  ice. 
This  energy  radiates  into  the  water,  forming  ambient  noise  when  the  events  are 
aggregated.  In  order  to  better  understand  the  Arctic  ambient  noise  generating 
mechanism,  two  levels  of  approaches  are  being  pursued  [9].  One  approach  is 
to  correlate  the  spectral  and  temporal  characteristics  of  ambient  noise  to  gross 
environmental  parameters  such  as  thermal  changes,  current  and  wind  stresses 
[33].  The  other  level  of  approach  is  to  look  into  the  individual  ice  cracking  events, 
and  treat  them  as  a  noise  generating  source  element,  which  forms  the  average 
ambient  noise  when  aggregated.  The  latter  approach  is  based  on  the  assumption 
that  the  central  Arctic  ambient  noise  is  dominated  by  the  radiation  due  to 
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mechanical  processes  in  the  ice  cover,  so  that  the  ambient  noise  is  aggregate  of 
individual  events  [9]  [34].  Therefore,  it  is  obvious  that  the  models  that  represent 
probable  source  mechanisms  of  a  single  event  need  to  be  developed  to  understand 
the  central  Arctic  ambient  noise  characteristics  as  an  aggregate  of  events. 

In  order  to  be  more  specific  about  the  issue  investigated  in  this  study  in  rela¬ 
tion  to  the  central  Arctic  ambient  noise,  the  processes  involved  in  the  generation 
of  Arctic  ambient  noise  can  be  categorized  as 

(1)  Development  of  environmental  stresses, 

(2)  Ice  plate  motion  induced  by  fracture  due  to  the  environmental  stresses, 

(3)  Radiation  from  the  ice  plate  into  water. 

The  first  phenomenon  has  been  extensively  studied  and  associated  with  the  ob¬ 
served  ambient  noise  [31]  [33].  The  radiation  mechanism  from  the  ice  plate  to 
water  is  also  studied  by  some  authors  [42]  [24]  [23]  [22].  However,  previous  pub¬ 
lications  are  mostly  speculative,  and  assume  certain  simple  source  types  in  the 
ice  plate.  Also,  it  is  not  clearly  understood  yet  how  the  environmental  stresses 
are  released  as  a  major  source  of  induced  ice  motion.  The  possible  candidates 
are  the  three  dominant  types  of  crack  [18]  .  Since  the  field  observation  of  such 
cracks  does  not  seem  to  be  feasible,  modeling  of  radiation  from  different  types 
of  crack  is  comparative  to  infer  the  source  mechanism  from  the  observed  signal 
of  the  ice  event. 

Another  motivation  of  the  modeling  will  be  the  expectation  that  an  un¬ 
derstanding  of  the  relationship  between  the  observed  acoustic  signal  and  the 
physical  process  of  sound  generation  will  lead  to  the  development  of  remote 
sensing  techniques  suitable  for  the  study  of  ice  behavior  and  properties  [10]. 

In  this  study,  the  scattering  is  not  included,  which  may  be  important  in  aver¬ 
age  ambient  noise  for  certain  frequencies  [22].  But,  the  observation  of  individual 
event  is  made  in  the  near  field,  where  presumably  the  scattered  field  is  much 
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sound  speed  in  the  water 


Figure  6.1:  Geometry  for  central  Arctic  environment. 


smaller  than  the  direct  field.  It  is  noted  that  the  attenuation  can  be  included 
by  using  complex  wave  number,  and  the  water  absorption  is  also  included. 

In  Section  6.1,  the  environmental  model  for  Central  Arctic  is  discussed.  In 
Section  6.2,  the  radiation  from  three  types  of  cracks  is  considered.  For  different 
types  of  sources,  the  radiation  pattern  and  temporal  characteristics  are  given 
and  discussed.  In  Section  6.4,  the  effect  of  the  non-compact  and  propagating 
crack  is  discussed.  Next,  in  Section  6.5,  it  is  stated  qualitatively  how  the  im¬ 
portant  parameters  affects  the  radiation  pattern  and  the  spectral  and  temporal 
characteristics.  In  Section  6.6,  the  effects  of  the  anisotropy  observed  in  the  sea 
ice  on  the  propagation  are  discussed.  Finally,  in  Section  6.7,  the  results  are 
summerized,  and  an  experiment  will  be  proposed  based  on  the  results. 
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6.1  The  Central  Arctic  Environment 


The  central  Arctic  environment  can  be  idealized  as  the  laterally  stratified  medium, 
assuming  smooth  boundaries1,  and  neglecting  the  scattering  due  to  ice  ridges 
and  roughness.  The  assumption  of  homogeneity  of  the  ice  can  be  justified  up 
to  the  frequency  of  wave  length  much  larger  than  the  grain  size  of  ice  separated 
by  brine  pocket.  However,  the  inhomogeneity  due  to  the  temperature  change 
along  thickness  is  expected  to  be  significant.  The  change  of  the  compressional 
and  shear  speed  due  to  the  seasonal  temperature  is  also  significant  [16],  but  the 
ice  is  assumed  to  be  homogeneous  for  the  purpose  of  studying  radiation  char¬ 
acteristics,  concentrating  on  the  source  mechanisms.  However,  it  is  noted  that 
there  is  no  limitation  in  obtaining  solutions  for  such  environmental  model. 

The  Arctic  environment  assumed  in  this  study  consists  of  6  layers  (Fig  6.1). 
The  first  layer  is  a  vacuum  half  space,  and  then,  a  homogeneous  elastic  ice 
layer  of  3  m  thickness.  The  next  three  layers  are  acoustic  layers  with  sound 
speed  gradient,  for  which  the  analytical  solution  can  be  obtained  in  the  form  of 
Airy  functions.  The  last  layer  is  a  fluid  half  space  with  constant  sound  speed. 
Although  the  bottom  can  be  replaced  by  a  sediment  layer  at  3000  m  deep,  the 
contribution  from  bottom  reflection  will  be  neglected. 


1The  scattering  due  to  rough  surface  has  been  treated  by  Kuperman  and  Schmidt  [21].  This 
scattering  effect  may  affect  significantly  the  average  ambient  noise,  but  it  is  neglected  in  studying 
the  direct  radiation  from  sources. 
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6.2  Radiation  from  Three  Modes  of  Compact 
Cracks 

In  this  section,  the  acoustic  emission  from  different  compact  sources  in  the  ice 
plate  will  be  studied  using  the  developed  numerical  code.  The  central  Arctic 
environment  in  Fig  6.1  is  considered,  and  a  compact  crack  is  assumed  to  be 
induced  at  1  m  from  the  top  of  the  ice.  The  source  time  series  and  spectrum 
for  the  seismic  moment  is  shown  in  Fig  6.2.  The  fault  orientation  parameters 
are  8  =  90°  and  <j>a  —  0°,  as  the  cases  considered  before.  Note  that  the  strike 
angle  <f>a  is  referred  to  x-axis,  i.e.  the  crack  surface  is  parallel  to  the  xz-plane  as 
shown  in  Fig  3.6.  The  following  cases  will  be  considered. 

Case  6.1  Compact  tensile  crack  with  receivers  at  z  —  60  m  :  The  results  are 
shown  in  Fig  6.3  for  radiation  pattern  and  synthetic  time  series.  The  receiver 
locations  are  same  as  for  the  previous  outputs  as  in  Fig  3.6. 

Case  6.2  Compact  tensile  crack  with  receivers  at  z  =  5  m  :  The  synthetic  time 
series  is  shown  in  Fig  6.4,  where  the  only  difference  from  the  Case  6.1  is  that 
the  receiver  depth  is  now  at  z  =  5  m  in  order  to  observe  the  radiation  from 
the  evanescent  first  antisymmetric  mode  in  the  floating  ice  plate.  The  inverse 
Hankel  integrands  for  the  different  receiver  depths  at  z  =  5  and  z  =  60  m  are 
given  in  Fig  6.6. 

Case  6.3  Compact  dip-slip  :  The  same  source  and  environment  as  in  the  Cases  6.1 
and  6.2  are  considered,  except  the  source  type  is  dip-slip.  The  radiation  pattern 
and  synthetic  time  series  are  given  in  Fig  6.7 

Case  6.4  Compact  strike-slip  :  The  same  source  and  environment  are  consid¬ 
ered,  except  the  source  type  is  strike-slip.  The  radiation  pattern  and  synthetic 
time  series  are  given  in  Fig  6.8 
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Figure  6.2:  Seismic  moment  (a)  source  function,  (b)  spectrum 
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Figure  6.3:  (a)  Radiation  pattern  for  /  =  100  Hz,  (b)  synthetic  time  series  for 
tensile  crack  with  dip  angle  6  =  90°  in  the  central  Arctic  environment. 
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Figure  6.4:  Synthetic  time  series  for  tensile  crack  with  dip  angle  6  =  90°,  re¬ 
ceivers  at  z  =  5  m,  r  =  300  m  at  every  30°  in  a  central  Arctic  environment. 

Two  kinds  of  graphical  display  will  be  used.  The  first  one  is  the  contour  plot 
of  the  pressure  field  at  60  m  depth,  in  Fig  6.3(a),  6.7(a),  and  6.8(a)  for  tensile 
crack,  dip-slip,  and  strike-slip,  respectively.  For  the  same  geometry  and  sources, 
the  synthetic  time  series  are  presented  in  Fig  6.3(b),  6.7(b),  and  6.8(b).  Since 
the  hydrophones  are  placed  in  the  water,  the  pressure  (negative  of  normal  stress) 
will  be  studied. 

Considering  the  tensile  crack  case,  the  pressure  field  in  the  water  at  60  m 
depth  shown  in  Fig  6.3  reflects  the  radiation  pattern  for  tensile  crack.  The 
time  series  in  Fig  6.3  (b)  shows  two  arrivals.  As  discussed  in  Stein  [42]  ,  the 
first  arrival  is  a  radiation  from  the  first  symmetric  mode  in  the  ice  plate,  often 
referred  as  the  longitudinal  mode.  The  second  arrival  is  the  so  called  acoustic 
mode,  i.e.  the  mainly  waterborn  arrival.  The  arrival  time  of  the  first  signal  can 
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DISPERSION  CURVES 


Figure  6.5:  Dispersion  curve  for  the  first  antisymmetric  mode  of  the  ice  plate, 
the  upper  line  is  for  the  group  velocity  and  the  lower  line  for  the  phase  speed. 

be  roughly  calculated  from  geometry  by  considering  the  critical  angle  at  the  ice- 
water  interface,  giving  0.121  sec.  The  second  arrival  time  is  0.212  sec.  Another 
mode  that  exists  is  the  contribution  of  the  radiation  from  the  first  antisymmetric 
mode  in  the  ice  plate,  which  is  often  referred  as  the  flexural  mode.  However, 
at  low  frequencies,  the  first  antisymmetric  mode  travel  subsonically  in  the  ice 
plate  with  an  evanescent  field  in  the  water.  Since  the  60  m  depth  is  2  wave 
lengths  deep  for  50  Hz,  which  is  the  lower  limit  of  the  frequency  components  of 
the  source  function  given  in  Fig  3.5,  this  mode  does  not  appear  in  the  synthetic 
time  series  shown  in  Fig  6.3(b).  In  order  to  see  the  dispersive  flexural  wave 
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Figure  6.6:  Modulus  of  inverse  Hankel  transform  integrand  for  Fourier  (a)  zeroth 
and  (b)  second,  cos  20,  orders  and  frequency  /  =  100  Hz,  at  z  =  5  m  and  z  = 
60  m. 
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Figure  6.7:  (a)  Radiation  pattern  for  /  =  100  Hz,  (b)  synthetic  time  series  for 
dip-slip  with  dip  angle  S  =  90°  in  a  central  Arctic  environment 
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radiation,  the  receiver  is  placed  at  z  —  5m,  for  which  the  synthetic  time  series 
for  pressure  is  shown  in  Fig  6.4.  At  this  point,  it  will  be  helpful  to  look  at  the 
modulus  of  the  integrand  of  inverse  Hankel  transform  for  each  Fourier  orders 
in  order  to  identify  how  the  energy  trapped  in  the  ice  plate  contribute  to  the 
pressure  below  the  ice.  The  moduli  of  integrand  for  depth  z  =  5  m  and  60  m 
axe  given  in  Fig  6.6.  The  modulus  of  the  integrand  for  Fourier  zeroth  order 
at  depth  z  =  60  m  has  two  peaks  at  horizontal  numbers  s  =  0.205  and  0.435. 
The  phase  velocity  of  the  first  peak  is  3065  m/sec,  and  represents  the  radiation 
from  the  first  symmetric  mode  in  the  ice  plate.  The  second  peak  appears  at  a 
speed  of  1444  m/sec,  and  corresponds  to  horizontal  acoustic  propagation  in  the 
water.  The  modulus  in  Fig  6.6  shows  another  peak  at  horizontal  wave  number 
s  =  0.675  corresponding  to  the  first  antisymmetric  mode  radiation.  This  mode 
has  a  phase  speed  of  931  m/sec.  However,  the  arrival  from  the  flexural  mode 
arrives  almost  same  time  as  the  acoustic  mode  in  Fig  6.4.  The  explanation  to 
this  arrival  time  difference  is  contributed  to  the  fact  that  the  group  velocity  of 
antisymmetric  mode  is  larger  than  its  phase  speed.  In  fact,  the  dispersion  curve 
of  the  first  antisymmetric  mode  of  ice  plate  can  be  determined  numerically  from 
the  Hankel  integrand  using  the  relation, 


ca  — 


du 
dk  ' 


(6.1) 


Eq  6.1  has  been  estimated  numerically  for  the  flexural  wave,  with  the  result 
shown  in  Fig  6.5  together  with  the  associated  phase  speed.  The  dispersion 
curve  shows  that  the  group  velocity  is  «  1450  m/sec.  and  the  phase  velocity 
is  ss  930  m/sec.  for  100  Hz,  in  agreement  with  the  observation  of  the  Hankel 
integrand.  In  the  high  frequency  limit,  the  phase  velocity  will  approach  to  the 
speed  of  the  interface  waves,  i.e.  Stoneley  waves  for  the  elastic-acoustic  and 
elastic-elastic  interfaces  or  Rayleigh  waves  for  the  elastic  half  space,  where  the 
wave  is  no  longer  dispersive.  Since  the  group  velocity  is  dispersive  for  the  present 
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cases,  the  waves  with  low  frequencies  arrive  late,  which  can  be  seen  in  Fig  6.4. 
For  the  second  Fourier  order  cos  29,  the  position  of  poles  are  the  same  but  the 
magnitude  is  about  half  of  the  zeroth  order,  which  explains  the  slowly  varying 
amplitude  along  the  azimuthal  angle  0  in  Fig  6.3  and  6.4.  Thus,  the  observation 
of  the  Hankel  integrand  gives  physical  insights  in  the  excitation  of  the  various 
modes,  and  their  corresponding  phase  speeds. 

For  the  dip-slip  crack,  the  pressure  field  and  time  series  are  given  in  Fig  6.7. 
The  amplitude  of  the  first  arrival,  i.e.  radiation  from  the  symmetric  mode  in  the 
ice  plate  is  of  much  smaller  amplitude  than  the  acoustic  arrival,  while,  for  strike- 
slip  (Fig  6.8),  the  amplitude  of  the  acoustic  mode  is  small.  This  is  contributed 
to  the  different  source  mechanisms.  For  dip-slip  with  dip  angle  8  —  90°  which  is 
presently  considered,  the  particle  motion  is  induced  in  vertical  direction  causing 
more  volume  fluctuation  at  the  ice-water  interface,  however,  for  strike-slip,  the 
particle  motion  is  in  horizontal  direction  exciting  the  compressional  wave  more 
effectively. 

Another  important  observation  can  be  made  from  the  transfer  functions  of 
each  crack  mode.  Denoting  the  source  function  in  Fig  6.2  as  s(f)  in  time  domain 
and  S(f)  in  frequency  domain  respectively,  and  the  system  transfer  function 
h(t )  and  H(f )  for  each  mode  of  cracks,  the  received  signal  at  a  receiver  position 
r  can  be  expressed  as 

r(t,r)  =  s(t)  *  h(t,  r) 

R(f,r)  =  ) 

in  the  time  and  frequency  domains.  The  transfer  functions  at  the  receiver  posi¬ 
tion  r  =300  m,  d  =  60  m,  and  6  =  30°,  are  shown  in  Fig  6.9.  The  vertical  axis 
represents 

Y  =  201og,0|B'(/)|.  (6.2) 
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The  transfer  function  characterizes  the  spectral  shape  of  the  received  signal  for 
an  impulse  source  function  of  a  particular  crack  type,  as  can  be  seen  in  Fig  6.9. 
However,  the  source  type  is  not  known  priori,  and  has  to  be  observed.  Thus,  the 
transfer  function  can  be  only  used  for  source  time  series  inversion  if  the  source 
type  is  known. 

It  has  been  shown,  from  the  figures  6.3,  6.7,  and  6.8,  that  the  radiation 
patterns  and  other  information  concerning  the  different  types  of  cracks  includ¬ 
ing  the  relative  strength  of  existing  modes,  can  be  obtained.  These  distinct 
characteristics  of  each  crack  mode  can  be  used  to  identify  the  dominant  source 
mechanisms  in  the  central  Arctic  environment  from  the  collected  data.  These 
discussions  are  found  in  the  following  section. 


6.3  Most  probable  source  types  and  correspond¬ 
ing  source  strength 

Most  Probable  Source  Types 

Among  the  three  examples  in  the  previous  section,  it  is  shown  that  the  radiation 
from  the  longitudinal  wave  in  the  ice  plate  is  much  stronger  than  the  acoustic 
mode  in  the  water  for  tensile  crack  and  strike-slip  with  dip  angle  6  =  90°,  which 
is  not  the  case  for  the  observed  signal  in  the  central  Arctic  environment2.  Thus, 
the  ratio  of  the  amplitude  of  the  radiation  from  the  longitudinal  wave  in  the  ice 
plate(Aj)  to  the  amplitude  of  the  acoustic  mode(Aa)  can  be  used  to  determine  the 
most  probable  source  types.  The  ratio,  are  shown  in  Fig  6.10  for  three  crack 
types  with  different  fault  orientations.  From  the  figure,  the  source  types  that 

2Discussion  with  Professor  Dyer 


118 


Dip  Angle 


Figure  6.10:  Ratio  of  the  amplitude  of  the  radiation  from  the  longitudinal  wave 
in  the  ice  plate(Aj)  to  the  amplitude  of  acoustic  mode(A0),  i.e.  for  three 
types  of  cracks  with  varying  dip  angle 

satisfy  the  relative  smallness  of  the  radiation  from  longitudinal  wave  compared 
to  the  acoustic  mode  are  dip-slip  with  dip  angle  S  fa  0°  and  90°  and  strike-slip 
with  dip  angle  6  fa  0°.  In  fact,  it  is  noted  that  the  mathematical  formulation  for 
strike-slip  with  dip  angle  6  =  0°  is  the  same  as  those  of  dip-slip  with  dip  angle 
<5  =  0°  and  90°. 

The  previous  study  by  Stein  [42]  used  acoustic  monopole  as  a  source  in  the 
ice  plate.  Due  to  the  nature  of  the  omnidirectional  volume  expansion,  the  ra¬ 
diation  from  the  longitudinal  wave  in  the  ice  plate  is  much  stronger  than  the 
acoustic  mode.  Study  by  Langley[24,23,22]  used  a  point  force[24,23]  and  tensile 
crack[22].  The  radiation  pattern  at  farfield  caused  by  tensile  crack  with  vertical 
fault  surface(dip  angle  6  =  90°)  shows  the  strong  radiation  from  the  longitudinal 
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wave  in  the  ice  plate  regardless  of  the  source  depth  at  frequencies  around  100 
Hz,  which  agrees  with  the  results  given  in  section  6.2. 


Source  Strength  Inversion 

Since  the  formulations  for  dip-slip  with  dip  angle  S  =  0°  and  90°,  and  strike-slip 
with  dip  angle  6  =  0°,  are  all  same,  dip-slip  with  dip  angle  90°  are  considered 
for  source  strength  inversion  from  the  pressure  in  the  water. 

The  pressure  amplitude  of  synthetic  time  series  in  Fig  6.11  at  range  r  = 
300  m  and  at  depth  z  =  60  m  in  the  horizontal  angle  6  =  30°  is  used  to  find  the 
source  strength.  The  seismic  moment  to  produce  pressure  1  Pa  at  the  specified 
receiver  position  is 

M  —  2.5  x  105iV  •  m 


Since  the  seismic  moment  is  M  =  Aun, 

,  2.5  x  105  K  , 

An  = - =  8.6  x  10~5  m3 


The  average  displacements  on  the  fault  surface  for  give  fault  surface  area  are 


A  m2 

u  mm 

10 

0.0086 

1 

0.086 

0.1 

0.86 

0.01 

8.6 

The  pressure  field  variation  depending  on  the  horizontal  angle  and  range  for 
200  Hz  is  shown  in  Fig  6.12. 


Although  the  amplitude  ratio  for  the  tensile  crack  with  dip  angle  6  =  0°, 
is  somewhat  large,  the  source  strength  will  be  given  since  the  amplitude  of 
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Figure  6.11:  Synthetic  time  series  in  the  horizontal  direction  30°  at  range 
r  =  300  m  and  at  depth  of  2  =  60  m  for  dip-slip  with  dip  angle  6  —  0° 
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Figure  6.12:  Field  pressure  variation  depending  on  the  horizontal  angle  and 
range  for  200  Hz.,  source  strength  is  4ir  N  •  m 


radiation  from  the  longitudinal  wave  is  affected  by  scattering  due  to  the  rough 
surface  of  ice  plate.  From  Fig  6.13,  the  seismic  moment  to  produce  pressure  1 
Pa  at  the  same  receiver  position  is 


M  =  7.5  x  106JV  •  m 


Since  the  seismic  moment  is  M  =  Au(X  +  2/z) , 


,  _  7.5  x  106  _4  , 

Au  —  — - - —  =  6.8  x  10  4  mr 

A  +  2  n 


The  average  displacements  on  the  fault  surface  for  give  fault  surface  area  are 
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Figure  6.13:  Synthetic  time  series  in  the  horizontal  direction  30°  at  range 
r  =  300  m  and  at  depth  of  z  =  60  m  for  tensile  crack  with  dip  angle  5  =  0° 
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6.4  Radiation  from  Propagating  Cracks 

The  effect  of  the  cracks  being  non-compact  is  important  for  the  frequencies  with 
wave  length  of  order  of  crack  length  or  less,  causing  frequency  dependent  di¬ 
rectivity  pattern.  In  reality,  the  crack  dimension  varies  from  the  order  of  grain 
boundary  length  to  order  of  103~4  m.  The  field  data  analysis  by  Farmer[lO] 
suggests  that  the  crack  could  be  non-compact,  and  the  estimated  crack  length 
range  between  27  and  63  m  for  maximum  length  based  on  the  analysis  of  four 
events.  For  the  crack  propagation  speed,  Mansinha[30]  showed  that  the  max¬ 
imum  crack  propagation  speed  for  a  medium  with  Poisson’s  ratio  of  0.25  is 
0.775c,  for  pure  shear  fracture  and  0.631c,  for  pure  tensile  fracture.  The  frac¬ 
ture  propagation  speed  increases  as  the  Poisson’s  ratio  increases.  Therefore,  for 
the  case  of  propagating  cracks  in  the  ice  plate,  the  propagation  speed  will  be 
assumed  to  1200  m/sec.  for  all  types  of  crack.  The  model  developed  in  Chapter 
4  is  applied  to  the  following  examples  with  a  set  of  typical  source  parameters 
to  investigate  the  effect  of  propagating  cracks.  The  radiation  from  propagating 
cracks  for  the  same  environment  (Fig  6.1)  is  considered.  The  source  types  are 
again  tensile,  dip-slip,  and  strike-slip  cracks  with  fault  dip  angle  6  =  90°  and 
strike  angle  <j>a  =  0°.  The  crack  geometry  is  considered  to  be  a  unilateral  line 
crack  with  length  10  m,  propagating  from  x  =  —  5m  to  +  5  m.  Since  the 
smallest  wave  length  A  =  cmin/fmax,  where  fmax  =  150  Hz  and  cmin  =  1800 
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m/ sec  for  shear  wave  speed,  is  greater  than  10  m,  the  spacing  between  source 
should  be  smaller  than  A/8,  which  corresponds  to  phase  difference  tt/4.  In  the 
following  examples,  the  distance  between  the  discretized  sources  are  taken  to 
be  1  m.  The  vertical  dimension  of  the  fault  surface  is  also  assumed  to  be  less 
than  A/8  to  be  treated  as  a  vertically  compact  line  crack.  In  order  to  consider 
a  crack  with  a  vertical  dimension  greater  than  A/8,  it  is  necessary  to  use  a  two 
dimensional  surface  crack  rather  a  line  crack.  Therefore,  the  propagating  cracks 
presently  considered  are  internal  cracks  centered  at  z  =  1  m. 

Case  6.5  Propagating  dip-slip  :  The  radiation  pattern  and  synthetic  time  series 
are  given  in  Fig  6.14 

The  pressure  field  for  dip-slip  for  100  Hz  harmonic  source  function  (Fig  6.14) 
is  compared  to  the  non-propagating  compact  crack  (Fig  6.7).  The  spatial  phase 
interference  causes  the  directivity  pattern  varying  rapidly  in  the  azimuthal  angle 
6.  The  dominant  radiation  direction  is  not  quite  clear  since  the  direction  of  major 
lobe  is  affected  not  only  by  the  radiation  pattern  of  the  particular  source  type, 
but  also  by  the  propagation  speed.  Since  the  crack  propagation  speed  is  subsonic, 
i.e.  less  than  the  sound  speed  in  the  water,  the  major  lobe  does  not  appear  in 
the  real  space.  However,  for  dip-slip,  the  acoustic  mode  is  dominant,  defining 
the  dominant  radiation  direction  more  clearly.  The  pressure  field  in  0°  and  180° 
directions  vanish  as  the  dip-slip  radiation  pattern  has  null  in  those  directions. 
This  is  an  indication  of  strong  dependence  of  the  radiation  of  propagating  crack 
on  the  source  types. 

An  important  observation  is  the  modulated,  spread  signal  in  the  time  se- 
ries(Fig  6.14(b)).  This  is  directly  contributed  to  the  modulation  introduced  by 
distributed  sources  over  the  fault  surface.  The  later  arrival  especially  in  the 
direction  of  0  =  150°  has  been  observed  in  the  crack  radiation  experiment  by 
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Figure  6.14:  (a)  Radiation  pattern  at  /  =  200  Hz,  (b)  synthetic  time  series  for 
dip-slip  with  dip  angle  6  =  90°  in  a  central  Arctic  environment. 
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Figure  6.15:  (a)  Synthetic  time  series,  and  (b)  spectrum,  in  the  horizontal  angle 
0  =  30°  for  moving  dip-slip 
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Figure  6.16:  (a)  Synthetic  time  series,  and  (b)  spectrum,  in  the  horizontal  angle 
9  =  90°  for  moving  dip-slip 
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Figure  6.17:  (a)  Synthetic  time  series,  and  (b)  spectrum,  in  the  horizontal  angle 
0  =  150°  for  moving  dip-slip 
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Savage  and  Mansinha  [37]  and  considered  as  the  stopping  phase. 


6.5  Parametric  Study  of  Crack  Radiation 

In  the  previous  sections  6.2  and  6.4,  the  radiation  pattern  and  temporal  char¬ 
acteristics  have  been  studied  for  a  set  of  source  and  environmental  parameters. 
Source  parameters  include  the  modes  of  crack,  fault  orientation  parameters, 
source  depth,  and  crack  propagation  speed  and  dimension  for  moving  cracks. 
The  effect  of  the  varying  source  parameters  on  the  radiated  field  will  be  impor¬ 
tant  to  identify  those  source  parameters,  since  the  interpretation  of  observation 
can  lead  us  to  better  understanding  of  the  source  mechanism. 

In  field  observation  of  the  acoustic  radiation,  the  signal  contains  the  infor¬ 
mation  on  the  source  mechanism.  The  information  consists  of  the  temporal  and 
spectral  characteristics  for  each  spatial  observation  positions.  The  complete  de¬ 
scription  of  the  field  due  to  all  source  parameters  can  take  up  the  rest  of  this 
thesis,  therefore,  only  the  qualitative  description  is  briefly  given.  It  is  considered 
to  be  another  research  area  to  be  studied  for  the  source  inversion  techniques. 

The  effect  of  crack  dimension  on  the  directivity  pattern  is  dependent  on  the 
frequency.  The  number  of  lobes  in  the  real  space  is  determined  by  the  relative 
length  of  the  crack  dimension  to  the  wave  length.  The  spectral  analysis  of  time 
series  at  spatially  distributed  sampling  positions  can  be  used  to  determine  the 
frequency  dependent  directivity  pattern  giving  the  estimation  of  crack  dimen¬ 
sion.  However,  this  problem  may  be  more  complicated  by  introducing  other 
parameters.  For  example,  the  crack  propagation  speed  might  vary  along  the 
crack  length,  and  the  crack  geometry  can  be  complex.  Although  the  theoretical 
synthesis  of  these  source  mechanisms  by  the  present  approach  is  possible,  the 
inversion  from  the  data  may  be  a  formidable  task,  unless  the  geometry  of  the 
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crack  is  known.  These  considerations  could  only  resolved  by  a  controlled  exper¬ 
iment,  where  the  measurement  of  source  parameters  is  possible.  The  design  of 
an  experiment  is,  therefore,  proposed  in  Section  6.7 

The  effect  of  source  depth  has  been  discussed  in  Langley  [23]  using  the  ver¬ 
tical  and  horizontal  forces  based  on  the  low  frequency  approximation.  It  shows 
that  the  importance  of  the  force  position  can  be  generally  neglected  for  low 
frequencies.  For  finite  frequencies,  however,  the  position  of  forces  is  important 
exciting  different  modes  of  waves  in  the  ice  plate.  For  example,  the  symmetric 
waves  are  likely  to  be  excited  when  the  source  depth  is  in  the  neutral  axis  of  the 
plate  including  the  fluid  loading.  Similarly,  the  antisymmetric  loading  excites 
the  antisymmetric  modes  more  effectively.  The  orientation  of  fault  surface  has 
significant  effect  on  the  radiated  field  for  the  tensile  and  shear  cracks  as  well  as 
the  point  forces.  The  results  for  the  varying  fault  orientation  parameters  are 
not  given,  but  their  effects  can  be  found  by  running  the  developed  numerical 
code.  In  fact,  the  inversion  of  the  fault  orientation  parameters  (refer  to  Fig  3.2) 
are  calculated  from  the  field  observation  using  various  inversion  techniques  in 
seismology. 

For  the  finite  frequencies,  where  the  thickness  of  the  ice  is  about  the  same 
order  or  greater  than  the  wave  lengths  of  the  frequencies  of  interest,  there  exist 
higher  modes  of  the  symmetric  and  antisymmetric  modes  in  the  ice  plate.  In 
turn,  each  modes  radiate  into  the  water  complicating  the  field.  In  higher  fre¬ 
quencies,  it  is  expected  that  the  interface  waves,  such  as  Rayleigh  and  Stonely 
waves,  propagate  along  the  interfaces  of  vacuum-ice  and  ice-water  interfaces. 
Since  these  waves  are  evanescent  spreading  cylindrically  (amplitude  ~  £),  the 
field  at  the  surfaces  of  ice  plate  is  dominated  by  them  in  the  long  range  propa¬ 
gation,  unless  there  is  scattering  due  to  the  rough  surface  of  the  ice  plate. 

Another  environmental  parameter  which  contains  useful  seismic  information 
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is  the  anisotropy  of  the  ice.  The  anisotropy  of  the  ice  is  caused  by  the  verti¬ 
cally  oriented  column  of  the  ice  separated  by  brine  pocket,  where  the  horizontal 
strength  of  the  ice  is  greater  than  the  vertical  direction.  The  distinguishing 
characteristics  of  the  transversely  isotropic  medium  is  the  separation  of  SH  and 
SV  waves  when  traveling  in  the  horizontal  direction.  This  is  discussed  in  the 
next  section. 


6.6  Anisotropy  of  Ice 

The  separation  of  SH  and  SV  waves  in  the  floating  ice  plate  in  the  central  Arctic 
ocean  has  been  reported  in  Hunkins  [16],  and  is  contributed  to  the  anisotropy 
of  the  ice.  The  anisotropy  is  caused  by  the  structure  of  the  ice,  where  the 
vertically  oriented  plates  separated  by  brine  brine  pockets.  Thus,  any  vertical 
shear  stress  would  act  along  the  many  lines  of  weakness  presented  by  the  vertical 
brine  pockets,  but  any  horizontal  shear  stress  would  meet  the  resistance  of  the 
interlocking  ,  randomly  oriented  grains.  Indeed,  the  measurement  of  static  ice 
strength  shows  that  the  lake  ice  as  well  as  the  sea  ice  are  not  isotropic  [16]. 

The  velocity  of  SV  and  SH  waves  vary  considerably  with  a  seasonal  change, 
which  is  largely  attributed  to  the  variation  of  ice  temperature  affecting,  again, 
the  elastic  constants.  Although  the  measurement  of  the  SH  wave  provides, 
when  observed  in  the  ice,  the  seismic  properties  of  ice  including  some  of  the 
elastic  constant,  this  wave  can  not  be  detected  in  the  water  because  the  SH 
wave  is  decoupled  from  the  vertical  displacement  of  ice,  so  that  the  water  wave 
is  not  excited.  Therefore,  it  is  concluded  that  the  SH  wave  does  not  radiate  any 
energy  into  the  water,  and  is  trapped  in  the  ice  plate  unless  there  is  scattering 
from  the  rough  surface. 

The  following  examples  for  the  isotropic  and  transversely  isotropic  media 
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will  show  that  the  effect  of  anisotropy  of  the  sea  ice.  The  compact  tensile  crack 
at  z  —  1  m  with  fault  orientation  parameters  6  =  90°  and  (f>t  =  0°  will  be 
considered. 

Case  6.6  Isotropic  ice  plate  in  the  central  Arctic  environment :  This  case  is  the 
same  as  the  Cases  6.1  and  6.2,  where  the  compressional  and  shear  velocities  are 
3500  m/sec  and  1800  m/sec,  respectively,  except  that  the  receivers  are  now  in 
the  ice  plate  instead  of  water.  Therefore,  the  field  parameters  are  calculated 
for  3  components  of  the  particle  velocities(Fig  6.18).  Also,  the  inverse  Hankel 
transform  integrands  for  each  components  of  field  parameters  and  Fourier  orders 
are  shown  in  Fig  6.21,  where  the  solid  line  is  for  the  isotropic  ice  plate  and  broken 
line  for  the  anisotropic  ice  plate  treated  in  the  next  Case  6.7. 

Case  6.7  Transversely  isotropic  ice  plate  in  the  central  Arctic  environment  : 
The  source  and  environmental  parameters  are  same  as  the  Case  6.6,  except 
that  the  ice  plate  is  now  transversely  isotropic.  The  transversely  isotropic  ice 
plate  consists  of  two  layers.  The  compressional  and  shear  velocities  are  4000 
m/sec  and  2300  m/sec  for  the  first  layer,  and  3000  m/sec  and  1300  m/sec  for 
the  second  layer,  respectively.  The  synthetic  time  series  for  the  3  components 
of  particle  velocities  are  given  in  Fig  6.19.  Also,  the  slowness  surfaces  for  the 
existing  waves  are  shown  in  Fig  6.20. 

First,  let  us  take  a  look  at  the  synthetic  time  series  for  Cases  6.6  and  6.7. 
For  the  vertical  displacements,  Case  6.6(Fig  6.18(a))  shows  that  the  earlier 
arrival  of  SV  Waves  of  flexural  waves  (second  arrival)  compared  to  that  of 
Case  6.7(Fig  6.19(a).  This  is  due  to  the  slow  SV  waves  in  the  transversely 
isotropic  medium  as  can  be  seen  in  Fig  6.21.  The  peaks  around  0.675  is  the 
contribution  from  the  SV  waves.  The  phase  speed  is 

CSV  =  ^  «  930  m.  (6.3) 
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Figure  6.18:  Synthetic  time  series  in  an  isotropic  ice  plate  in  the  central  Arctic 
environment,  (a)  vertical  displacement,  (b)  x-direction  displacement. 
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Figure  6.19:  Synthetic  time  series  in  a  transversely  isotropic  ice  plate  in  the  cen¬ 
tral  Arctic  environment,  (a)  vertical  displacement,  (b)  x-direction  displacement. 
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Norm,  horizontal  slowness 


Figure  6.20:  Slowness  surface  of  P(solid  line),  SV (dotted  line),  and  S'# (dashed 
line)  waves. 

However,  it  is  noted  that  the  group  velocity  is  much  faster  than  the  phase  speed 
for  the  SV  wave  giving  «  1450  m/sec  from  Fig  6.5,  which  explains  the  arrival 
time  of  approximately  0.212.  The  second  peak  of  broken  line  in  Fig  6.21(a) 
for  the  transversely  isotropic  ice  plate  represents  SV  wave  in  the  transversely 
isotropic  ice  plate  traveling  slower  than  the  SV  wave  in  the  isotropic  ice  plate. 
The  second  arrival  in  the  synthetic  time  series  in  Fig  6.18(b)  is  the  SH  wave 
in  the  isotropic  ice  plate,  which  travels  slightly  slower  than  the  SH  wave  in  the 
transversely  isotropic  ice  plate.  This  can  be  explained  again  from  the  position 
of  the  peaks  for  both  media  in  Fig  6.21(d)  and  (f).  The  second  peak  with  solid 
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Figure  6.21:  Modulus  of  Hankel  transform  integrand  for  isotropic  medium  (solid 
line),  and  transversely  isotropic  medium  (dotted  line),  (a)  and  (b)  show  the 
vertical  displacements  for  Fourier  orders  zero-th  and  cos  20  orders,  respectively, 
(c)  and  (d)  for  the  x-direction  velocity,  and  (e)  and  (f)  for  the  y-direction  velocity. 
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line  represents  the  wave  number  for  the  isotropic  ice  plate  and  that  with  broken 
line  for  the  transversely  isotropic  ice  plate,  explaining  eaxlier  arrival  of  the  SH 
wave  in  the  transversely  isotropic  ice  plate  than  that  of  the  isotropic  ice  plate. 
The  phase  speed  of  the  SH  waves  are 

cSh  =  j  «  1800  m,  (6.4) 

which  is  almost  the  same  as  the  medium  SH  wave  speed  of  the  ice  (Fig  6.20). 
This  is  because  the  SH  wave  is  completely  decoupled  from  the  waves  existing 
in  the  water.  The  exact  values  of  the  phase  speed  of  .S' if  wave  can  be  found  in 
Fig  6.20  for  the  different  directions  of  propagation. 

The  arrival  of  the  compressional  waves  for  both  media(Fig  6.18  and  6.19)  is 
not  much  affected  as  can  be  seen  in  Fig  6.21  for  the  modulus  of  Hankel  integrand 
at  the  top  of  the  ice. 


6.7  Proposed  Crack  Radiation  Experiment 

In  spite  of  the  large  amount  of  ambient  noise  data  collected  over  the  years,  no 
data  sets  available  are  directly  suited  for  studying  the  individual  events  to  con¬ 
firm  existing  crack  radiation  hypotheses.  This  is  due  to  several  factors  including 
the  limited  spatial  sampling  imposed  by  economy  and  the  fact  that  most  exper¬ 
iments  have  focused  on  the  low  frequency,  overall  ambient  noise,  requiring  large 
apertures.  Also,  the  layout  of  sensors  requires  a  certain  arrangement  suited  for 
studying  the  radiation  patterns  from  the  sources.  Therefore,  in  this  section,  a 
series  of  experiment,  i.e.  laboratory  and  field  experiments,  will  be  proposed. 

The  objective  of  the  proposed  laboratory  and  field  experiments  are 

1.  Laboratory  Experiment 

•  Validation  of  model  with  the  controlled  sources. 
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•  Development  of  inverse  technique. 


2.  Field  Experiment 

•  Source  parameters  inversion  from  the  observed  data  including  source 
location,  crack  type,  crack  dimension,  and  propagation  speed. 

The  data  obtained  in  the  laboratory  will  be  analyzed  to  show  the  radiation 
pattern  in  the  ice  and  water,  and  used  to  test  the  existing  analytical  and  nu¬ 
merical  models.  This  will  be  done  by  investigating  the  directivity  pattern  as 
well  as  the  temporal  characteristics  for  different  modes  of  cracking  or  different 
forcing  mechanisms.  Also,  the  effects  on  the  pressure  field  in  the  water  of  the 
non-compactness  of  a  source  with  finite  crack  propagation  speed  in  the  floating 
ice  sheet  will  be  studied.  It  is  the  expected  outcome  of  this  proposed  experiment 
to  develop  the  inversion  techniques  from  the  known  source  parameters  and  the 
observed  data.  Once  the  spectral  and  temporal  characteristics  of  each  fracture 
modes  are  identified  by  fitting  data  through  our  analytical  models,  these  models 
can  be  applied  to  the  data  collected  in  the  central  Arctic  environment  in  various 
field  experiments.  For  example,  radiation  patterns  for  each  mode  of  cracking 
can  be  used  to  identify  which  modes  of  cracking  are  dominant,  and  how  they 
radiate  energy  into  water  in  the  central  Arctic  environment.  Thus,  better  un¬ 
derstanding  of  the  individual  crack  radiation  events  can  lead  us  to  refined  Arctic 
Ocean  ambient  noise  model. 

The  model  developed  in  this  research  served  as  a  guide  to  design  the  ex¬ 
periment  and  gives  ideas  about  what  is  to  be  sought  from  the  analysis  of  data. 
The  idea  is  to  relate  the  characteristics  of  the  observed  data  to  the  features 
predicted  by  the  model,  yielding  the  source  parameters  of  the  event.  The  source 
parameters  thus  obtained  are  directly  related  to  the  physical  source  mechanisms. 
Following  are  discussions  about  what  aspects  are  to  be  sought  from  the  data. 
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The  details  of  these  aspects  are  found  in  the  sections  6.2,  6.3,  and  6.4. 

In  the  frequency  domain,  the  spectral  level  at  each  spatially  distributed  sen¬ 
sor  is  plotted  to  display  the  directivity  pattern  at  each  frequency.  This  will  be 
done  for  a  known  mode  of  cracking  with  crack  parameters  such  as  dimension, 
direction  and  propagation  speed.  This  directivity  can  be  used  to  identify  the 
source  type  and  crack  size  effect.  This  methodology  has  been  developed  and 
widely  used  in  seismology.  Also  the  slope  of  roll-off  in  spectrum  for  each  sensor 
location  can  be  used  for  identifying  the  source  time  series.  The  temporal  char¬ 
acteristics,  including  the  duration  of  observed  signal  and  the  stopping  phase, 
can  be  related  to  the  dimension  of  the  crack. 

6.7.1  Laboratory  Experiment 

In  order  to  investigate  the  crack  radiation  problem,  it  is  necessary  to  minimize 
scattering  due  to  rough  surface,  pre-existing  cracks,  or  ridges  so  as  to  concen¬ 
trate  on  the  source  inverse  problem.  In  the  Arctic  Ocean,  scattering  due  to 
inhomogeneity  of  the  ice  such  as  rough  surfaces,  pre-existing  cracks  and  ridges 
are  complicating  factors  obscuring  the  crack  radiation  problem.  It  is  therefore 
essential  initially  to  perform  a  laboratory  experiment  which  eliminates  these 
complicating  factors. 

The  laboratory  experiment  should  include  direct  measurement  of  crack  pa¬ 
rameters,  such  as  size,  propagation  speed  and  direction.  Crack  size  and  direction 
can  be  observed  after  the  event,  however  the  propagation  speed  is  hard  to  mea¬ 
sure  because  of  its  highly  transient  nature.  Here,  a  high  speed  camera  could  be 
used  provided  that  the  precise  position  and  time  of  crack  formation  is  known. 
Another  way  is  to  measure  the  time  history  of  the  displacement  along  and  across 
the  crack,  e.g.  by  means  of  accelerometers  or  strain  gauges. 

Inducing  a  desired  mode  and  size  of  crack  without  interfering  with  the  wave 
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field  is  difficult,  and  is  another  research  area  being  studied  3.  Although  the 
mechanical  properties  of  ice  are  rather  complex,  it  is  known  that  the  ice  is  much 
stronger  in  compression  than  in  tension.  It  is  possible  that  this  fact  may  be 
exploited.  In  order  to  ensure  the  crack  is  formed  at  the  desired  position  with 
certain  direction,  the  prescribed  crack  surface  is  cut  so  as  to  reduce  the  cross 
sectional  area.  These  aspects  of  inducing  different  modes  of  cracks  can  be  found 
in  the  ice  mechanics  literature  [8]. 

The  field  parameters  to  be  observed  include  the  pressure  in  the  water  as  well 
as  the  displacement  in  the  ice  plate.  Since  the  ice  is  an  elastic  plate,  supporting  a 
large  number  of  Lamb  modes,  it  is  expected  that  valuable  information  concerning 
the  cracking  process  can  be  obtained  by  measuring  the  ice  motion  in  addition 
to  the  sound  field  in  the  water.  The  other  observables  are  the  crack  parameters, 
such  as  the  volume  expansion  for  tensile  cracks,  the  seismic  moment  for  shear 
cracks,  the  dimension  of  the  crack,  and  the  crack  propagation  speed.  These  crack 
parameters  can  be  best  observed  by  measuring  the  particle  motion  at  the  vicinity 
of  the  crack  surface.  In  order  to  measure  the  time  history  of  displacement  at  the 
crack  surface,  it  is  suggested  to  use  high  speed  camera  with  a  resolution  enough 
to  observe  the  moving  crack  tip.  Another  way  is  to  place  the  accelerometers  or 
strain  gauges  along  the  prescribed  crack  allowing  for  localization  of  the  moving 
crack  tip. 

Temporal  Sampling 

The  sampling  frequency  of  the  sensors  is  determined  based  on  the  size  of  the 
crack  to  be  induced,  since  we  are  looking  at  the  radiation  pattern  of  non-compact 
source.  Non-compactness  of  the  crack  is  determined  by  the  size  of  the  crack 
relative  to  the  wave  length.  The  non-compact  source  has  phase  interference  due 

3  Discussion  with  Professor  T.  Wierzbicki 
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to  its  finite  size  of  source  distribution.  This  spectral  behavior  of  the  pressure 
field  (  or,  stress  field  in  the  ice  )  show  the  periodic  nulls  and  peaks  depending 
on  the  frequency,  which  is  called  frequency  dependent  directivity.  In  order  to 
identify  this  behavior,  the  observation  of  data  with  high  Nyquist  frequency  is 
necessary.  For  this  purpose,  the  sampling  frequency  of  10  KHz  is  suggested. 
Above  approximately  10  KHz,  the  assumption  of  homogeneity  of  ice  does  not 
hold,  since  the  roughness  and  grain  sizes  axe  greater  than  the  wave  length. 

Spatial  sampling 

Although  the  presently  developed  solution  technique  is  not  restricted  to  the 
far  field,  the  far  field  measurement  have  some  advantages  in  applying  to  the 
average  ambient  noise,  so  that  it  is  recommended  to  place  the  receivers  and 
accelerometers  in  the  far  field.  If  the  sensors  are  in  the  near  field,  the  radiation 
pattern  will  be  range  dependent  so  that  the  directivity  pattern  is  hard  to  identify. 
Moreover,  most  analytical  solutions  are  valid  only  in  the  far  field.  Conditions  for 
the  far  field  approximation  are  /  >C  r  for  amplitude,  and  ^  f  for  phase,  where 
r  is  the  distance  from  crack  to  observation  point.  These  conditions  reduces  to 

l  *C  r  (6.5) 


2n 


(6.6) 


where  c  is  the  medium  phase  speed,  /  is  the  frequency  of  interest,  k  is  wave 
number,  l  is  the  crack  dimension,  and  r  is  the  receiver  range. 

Based  on  the  considerations  presented  in  the  previous  sections,  it  is  recom¬ 
mended  that  the  spatial  sampling  is  at  every  30°  at  least  in  the  horizontal  and 


vertical  directions  to  identify  the  spatial  directivity  due  to  the  compact  direc¬ 
tional  sources,  as  well  as  the  frequency  dependent  directivity  due  to  the  moving 
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cracks. 


Design  Of  Array  Of  Sensors 

Based  on  the  considerations  presented  in  the  previous  sections,  an  array  design 
will  be  proposed.  To  measure  the  displacement  field  in  the  ice,  16  tri-axial 
accelerometers  need  to  be  placed  in  the  ice(Fig  6.22).  The  distance  from  the 
crack  to  the  accelerometers  r  is  recommended  to  be  as  large  as  the  dimension 
of  the  laboratory  allows,  so  that  the  longitudinal  and  flexural  modes  in  the  ice 
plate  can  be  observed  separately.  When  the  accelerometers  are  placed  at  r  « 
150  m,  the  10  Hz  frequency  components  will  be  separated  by  a  time  delay  of  one 
wave  length  traveling  time,  i.e.  0.1  seconds  (refer  to  Fig  6.5),  for  flexural  mode. 
Similarly,  a  time  delay  of  10  wave  lengths  of  traveling  time  will  be  achieved  for 
the  100  Hz  frequency  component.  If  the  distance,  r,  is  15  m,  the  separation 
time  will  be  reduced  by  a  factor  of  10.  When  the  separation  of  existing  modes 
is  critical,  the  scaled  experiment  with  artificial  is  suggested  to  reduce  the  size  of 
laboratory.  Alternatively,  the  lake  ice  can  be  used,  where  the  sensors  are  placed 
far  enough  to  separate  the  different  modes. 

To  measure  the  pressure  field  in  the  water,  20  hydrophones  are  to  be  sus¬ 
pended  in  the  water(Fig  6.23).  The  hydrophones  can  be  placed  about  same 
distance  as  the  accelerometers  in  order  to  separate  the  longitudinal  mode  from 
the  acoustic  mode.  The  depth  of  hydrophones,  however,  need  to  be  greater  than 
half  the  wave  length  at  the  frequency  of  interest,  so  that  the  flexural  wave  does 
not  interfere  with  the  other  modes.  For  example,  when  the  frequency  of  interest 
is  100  Hz,  the  depth  of  hydrophones  should  be  greater  than  d  «  3.5  m. 

The  crack  parameters  are  found  from  the  time  history  of  the  particle  motion 
at  the  crack  surface.  To  measure  the  displacement  along  and  cross  the  crack,  tri- 
axial  accelerometers  are  placed  as  indicated  in  Fig  6.24.  As  mentioned  earlier, 
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Figure  6.23:  (a)  Top  view  of  hydrophone  array  and  (b)  Cross  section  along 
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Figure  6.24:  Tri-axial  accelerometers  placed  along  the  prescribed  crack  surface 
to  measure  the  displacement  along  and  across  the  crack  surface. 

a  high  speed  camera  is  recommended  along  with  the  accelerometers. 

6.7.2  Field  Experiment 

Once  the  controlled  experiment  in  the  laboratory  is  carried  out,  and  the  features 
of  the  data  are  properly  explained  and  matched  with  the  model  developed,  a 
field  experiment  should  be  performed.  The  difference  between  the  laboratory 
and  field  experiment  will  be  the  randomness  of  the  event.  The  likelyhood  of 
observing  events  within  the  circular  array  of  radius  300  m,  in  which  we  can 
obtain  data  points  in  all  directions,  can  be  estimated  from  previous  experiments. 
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For  example,  in  the  FRAM  IV  study  [44],  it  has  been  shown  that  for  the  center 
annulus  of  30  square  Km,  there  were  91  events  per  30  square  Km  per  662  minutes 
of  observation,  or  approximately  0.3  events  per  square  Km  per  hour.  Since  our 
circular  array  has  0.28  square  Km,  we  expect  about  0.1  events  per  hour  within 
the  circular  array.  Therefore,  the  expected  duration  of  the  experiment  should  be 
10  Hrs  to  observe  an  event  within  the  circle.  However,  since  the  signal  to  noise 
ratio  of  an  event  depends  on  the  distance  from  the  event  to  the  hydrophones, 
we  might  observe  more  events  in  the  small  circular  area  enclosed  by  an  array  of 
hydrophones. 


6.8  Conclusion 

The  developed  model  has  been  applied  to  the  various  examples  to  demonstrate 
the  solution  technique,  and  to  find  a  number  observations  for  the  radiation 
pattern  and  temporal  characteristics  depending  on  the  source  and  environmental 
parameters.  The  source  parameters  can  be  the  types  of  source,  or  mode  of  crack, 
dimensions  of  crack  as  well  as  geometry  (fault  surface  orientation  parameters 
A ,<f>s,6,  refer  to  Fig  3.2),  propagation  speed  and,  source  depth.  Among  the 
environmental  parameters,  the  parameters  of  importance  are  thickness  of  the 
plate,  the  mechanical  properties  of  the  medium,  such  as  the  anisotropy  and 
inhomogeneity.  The  change  of  the  radiation  pattern,  and  temporal  and  spectral 
characteristics  are  discussed. 

First,  the  radiated  field  from  the  three  types  of  fracture  show  the  distinct  di¬ 
rectivity  pattern  when  observed  in  the  water,  depending  on  the  fault  orientation 
parameters.  The  presence  of  the  ice-water  interface,  coupling  the  P  and  SV 
waves,  complicates  the  radiation  pattern  as  well  as  the  time  series  separating 
the  contributions  from  the  radiation  of  different  modes  existing  in  the  ice  plate. 
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The  relative  strength  between  the  radiation  from  the  acoustic,  antisymmetric, 
and  symmetric  modes  in  the  ice  plate  for  the  three  types  of  fracture  is  used  to 
infer  the  most  probable  source  types  and  corresponding  dip  angle  based  on  the 
observation  that  the  acoustic  mode  is  much  stronger  than  the  radiation  from  the 
longitudinal  wave  in  the  ice  plate.  It  is  found  that  the  shear  cracks,  i.e.  dip-slip 
with  dip  angle  6  «  0°  and  90°  and  strike-slip  with  dip  angle  6  «  0°.,  are  the 
most  probable  source  types.  The  corresponding  source  strength  is  related  to  the 
field  pressure  in  terms  of  the  seismic  moment  and  the  average  displacement  over 
specified  fault  surface  area. 

The  temporal  characteristics  of  the  non-compact  and  propagating  crack  can 
be  characterized  as  the  modulated  signal,  of  which  cross-correlation  between 
channels  (receiver  positions)  are  very  poor.  Another  noticeable  result  is  the 
rapid  change  of  directivity  in  the  varying  angle.  The  larger  the  dimension  of  the 
crack  is,  the  more  peaks  and  nulls  of  the  radiation  pattern  appears  to  give  rapid 
variation  with  angle. 

Next,  the  effect  of  anisotropy  of  the  sea  ice  is  discussed.  The  sea  ice  is  shown 
to  be  anisotropic,  where  the  SH  and  SV  waves  separates.  The  additional  mode 
existing  due  to  anisotropy,  however,  does  not  affect  the  radiation  pattern,  since 
the  SH  wave  in  the  ice  plate  does  not  excite  the  vertical  displacement  at  the 
ice-water  interface  so  that  the  water  wave  is  decoupled  from  the  SH  wave. 

The  sources  in  the  ice  plate  as  an  element  of  Arctic  ambient  noise  can  con¬ 
tribute  to  the  average  noise  either  as  the  direct  radiation  or  as  the  scattered  field 
due  to  the  inhomogeneity  and  rough  boundary  of  the  ice.  In  certain  frequencies, 
the  different  types  of  crack  modes,  or  physical  mechanism  can  be  important. 
Further,  the  ambient  noise  forming  process  in  other  frequency  regimes  may  be 
completely  different,  as  well  as  the  responsible  environmental  stresses.  The 
application  of  this  model,  therefore,  should  be  made  with  the  careful  physical 
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insights  to  the  overall  ambient  noise  frequencies  as  well  as  the  corresponding 
generation  mechanisms  and  driving  forces. 


Chapter  7 
Conclusions 


7.1  Summary 

The  summary  of  this  thesis  emphasizing  the  original  contributions  can  be  made 
under  two  categories.  One  is  the  development  of  the  efficient  numerical  algo¬ 
rithm  as  well  as  the  corresponding  solution  techniques  for  radiation  from  the 
compact  and  non-compact  directional  seismic  sources  in  a  laterally  stratified 
medium,  combining  the  existing  numerical  code  (SAFARI)  with  source  repre¬ 
sentation  in  an  unbounded  medium.  The  compact  source  representation  by 
Keilis-Borok  (1950),  reviewed  in  Sato  [36],  is  treated  in  Chapters  2,  and  serves, 
in  nature,  as  tutorial  and  background  for  the  further  development  in  the  Chap¬ 
ter  3.  In  Chapter  3,  which  is  the  kernel  of  this  thesis,  the  representation  of 
seismic  sources  is  extended  to  the  formulation  of  tensile  crack  in  addition  to 
the  shear  fault  formulation,  and  transformed  to  the  cylindrical  coordinates  and 
arranged  in  such  a  way  that  the  formulation  can  be  applicable  to  the  global  ma¬ 
trix  method.  Consequently,  the  source  terms  are  incorporated  into  the  global 
matrix  method,  treating  the  multiple  azimuthal  Fourier  orders  simultaneously. 
In  Chapter  4,  the  superposition  model  to  treat  the  radiation  from  non-compact 
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array  and  propagating  crack  model  taking  advantage  of  numerical  method,  also, 
has  been  developed  eliminating  the  numerical  inefficiency  and  convergence  prob¬ 
lems  caused  by  the  expansion  of  the  field  with  Bessel  function[38].  This  non¬ 
compact  model  is  of  great  practical  importance  with  potential  applications  to 
the  radiation  from  an  array  of  airgun  or  explosions  for  seismic  exploration  and 
propagation.  Another  problem  treated  in  Chapter  5  is  the  propagation  effect 
of  the  transversely  isotropic  medium,  of  which  applications  are  found  in  the  sea 
and  lake  ice,  and  the  periodic  finely  layered  sediment.  Since  the  theory  for  the 
transversely  isotropic  medium  is  well  established,  the  focus  is  on  the  formulation 
of  the  compatible  boundary  conditions  with  the  global  matrix  approach. 

The  other  category  is  the  outcome  of  the  application  of  the  developed  model 
to  the  physical  system.  The  developed  model  is  basically  applied  to  the  canoni¬ 
cal  problems  and  to  the  environmental  model  of  the  central  Arctic.  In  Chapter 
3,  the  canonical  examples  of  radiation  from  the  three  different  modes  of  crack  in 
an  homogeneous  unbounded  medium  are  discussed.  The  spatial  and  temporal 
characteristics  of  the  radiated  field  are  distinctively  different  for  the  three  modes 
of  crack,  i.e.  tensile  crack,  dip-slip,  and  strike-slip.  The  detailed  discussions  are 
given  in  Chapter  3.5.  The  application  of  non-compact  and  propagating  crack 
model  is  found  in  Chapter  4,  where  the  radiated  field  are  characterized  by  the  fre¬ 
quency  dependent  directivity  pattern,  and  the  temporal  characteristics  include 
the  existence  of  the  stopping  phase,  the  modulation  of  signal  by  the  dimension 
of  the  crack  surface,  and  poor  cross  correlation  due  to  the  phase  interference 
due  to  the  crack  dimension.  The  results  of  the  application  to  the  central  Arctic 
environment  shows  more  complications  introduced  by  a  wave  guide  in  addition 
to  the  general  characteristics  mentioned  above  for  the  canonical  problems.  The 
coupling  of  the  P  and  SV  waves  at  the  ice-water  interface  introduces  the  dis¬ 
persion  of  each  components  of  waves,  and  the  corresponding  radiated  field  in 
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the  water.  For  the  directional  sources,  the  SH  wave  is  excited  and  trapped  in 
the  ice  plate,  since  the  SH  wave  decouples  from  the  P  and  SV  waves  and  there 
exists  only  the  P  wave  in  the  water.  When  the  rough  boundary  with  scattering 
effect  is  considered,  the  contribution  of  the  scattered  field  from  SH  and  SV 
to  the  average  ambient  noise  could  be  significant,  of  which  effect  needs  to  be 
studied  in  a  organized  way.  Finally,  the  effect  of  anisotropy  of  the  ice  reported 
by  Hunkins[l6]  is  simulated  to  give  the  separation  of  the  SV  and  SH  waves, 
which  contains  the  information  on  the  elastic  constants  in  different  directions  as 
well  as  the  dispersion  relation  for  the  given  geometry. 


7.2  Discussions  and  Suggestions  for  Future  Stud¬ 
ies 

Solution  Techniques 

The  applicability  of  presently  developed  solution  technique  is  limited  to  the 
range  independent  medium,  which  is  often  referred  as  the  laterally  stratified 
medium,  the  vertically  varying  medium,  or  the  laterally  homogeneous  medium. 
This  environmental  model  is  of  the  great  practical  importance  for  many  physical 
problems,  such  as  the  seismic  propagation  and  underwater  acoustic  propagation, 
where  the  medium  can  be  treated  as  laterally  homogeneous.  As  the  more  accu¬ 
rate  modeling  becomes  necessary,  the  laterally  inhomogeneous  medium  needs  to 
be  modeled,  especially  for  high  frequencies  where  the  effect  of  the  local  inhomo¬ 
geneity  is  significant,  of  which  example  can  be  found  in  the  scattering  from  the 
ridges  in  the  central  Arctic  environment  affecting  the  directly  observed  signal  as 
well  as  the  average  ambient  noise  as  result  of  long  range  scattering  effect. 
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Applications 


Application  of  the  developed  code  can  be  found  in  seismology  as  well  as  in 
underwater  acoustics.  The  applications  include  the  seismic  exploration,  the  field 
solution  to  the  radiation  from  the  sonar,  and  the  design  of  the  efficient  volume 
transducer  array  with  focusing  and  steering.  For  given  environment,  the  source 
discrimination  can  be  made  from  the  spectral  and  temporal  characteristics  of 
the  observed  data. 

As  a  source  inverse  problem,  the  developed  code  can  be  instrumental  in 
studying  the  source  mechanisms  of  fracture  in  the  floating  ice  plate  in  the  central 
Arctic  environment  as  an  element  of  the  average  ambient  noise.  The  directivity 
pattern  of  the  compact  and  non-compact  fracture  and  the  temporal  character¬ 
istics  including  the  dispersion  and  separation  of  different  modes  existing  in  the 
ice  plate  as  well  as  the  field  in  the  water  for  different  types  of  fracture  have  been 
found  to  be  readily  applicable  to  a  proper  set  of  data.  For  continuing  research 
effort,  controlled  crack  radiation  experiment  is  proposed  in  Chapter  6.6  and 
given  with  some  insights  to  the  layout  and  the  measurement  of  field  parameters 
based  on  the  preliminary  outcome  of  the  model.  The  experiment  will  be  used 
in  identifying  the  dominant  crack  mechanisms,  the  radiation  and  propagation 
effects  that  might  be  important  for  different  environment.  The  broader  issues, 
to  which  the  developed  model  and  crack  radiation  experiment  is  instrumental 
in  relation  to  the  central  Arctic  ambient  noise,  are 

•  in  what  frequency  band  is  the  crack  radiation  most  important  ? 

•  what  other  mechanisms  are  dominant  in  other  frequencies  ? 

Regarding  to  the  first  question,  only  the  field  observation  can  be  source  of  in¬ 
formation  that  leads  to  the  answer.  Again,  the  problem  is  complicated  by  the 
varying  environment,  which  modifies  the  natural  frequencies  of  the  environ¬ 
mental  system.  However,  the  properly  scaled  and  conditioned  experiment  and 
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analysis  is  the  most  probable  candidate  as  means  of  the  answering  the  raised 
questions. 

Therefore,  it  is  concluded  that  two  lines  of  research  be  pursued  towards 
better  understanding  of  the  crack  radiation.  First  relates  the  proper  modeling 
of  the  physical  system,  as  being  undertaken  by  many  researchers.  These  include 
the  propagation  effect  in  the  range  dependent  environment  including  scattering, 
and  other  kinds  of  sources,  such  as  gravity  waves  at  extreme  low  frequency, 
radiation  and  propagation  at  high  frequency  where  the  medium  can  no  longer 
be  treated  as  homogeneous,  and  so  on.  Second,  the  experimental  approach  using 
the  field  and  laboratory  data  should  be  used  to  define  the  dominant  mechanisms 
as  well  as  confirm  the  prediction  made  by  modeling,  which  will  eventually  lead 
to  better  understanding  of  the  physical  mechanism  and  to  constructing  a  model 
to  predict  the  necessary  informations. 
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Appendix  A 

Potentials,  Displacements  and 
Stresses  of  Various  Sources  with 
an  Arbitrary  Direction  in 
Cylindrical  Coordinates 


The  potentials,  displacements  and  stresses  caused  by  various  sources  are  pre¬ 
sented  in  this  section.  The  potentials  are  related  to  displacements  in  the  follow¬ 
ing  manners. 

u  =  V^  +  VxVx  (A.l) 

=  V<£  +  V  x  V  x  (0,0,  A)  +  V  x  (0,0,^)  (A.2) 

The  potentials  are  summed  over  each  Fourier  order  m,  where  the  m-th  order  are 
denoted  by  superscript  m.  Using  the  potentials  defined  in  Eq  A.l,  the  potentials 
are 

OO 

M,z)  =  'E<T(r,z) 

m= 0 


cos  mO 
sin  mO 
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(A.3) 


oo 

cos  mO 

A(r,»,*)  =  E‘*(v) 

m—0 

sin  mO 

• 

oo 

sinm# 

=  I 

m= 0 

—  cos  md 

Using  the  definition  in  Eq  A.2,  the  potentials  are 


OO 

cosm0 

4>{r,0,z)  =  ^4>m{r,z) 

m=  0 

sinm0 

oo 

cos  m# 

VvM.z)  = 

m—0 

sin  mO 

oo 

sin  mO 

ijje(r,0,z)  =  Y^i>T{r,z) 

m—0 

—  cos  m.0 

oo 

cos  mO 

«r,M  =  £  C(r.») 

m=0 

sinra# 

(A.4) 


A.l  Point  Force  in  an  Arbitrary  Direction 

Denoting  the  force  vector  F  =  Fxez + Fvey + Fzez  as  a  source  term,  the  potentials 
and  the  field  parameters  are  as  follows. 


Potentials 


<i>  = 


exwt 

4wpu)2 


Fz$aJ0(sr) 
+Fxcos$  sJi(sr) 
+FvsinO  sJi(sr) 
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A  = 


0  = 


eiwt  r° 
4npu>2  Jo 

eiut  n 
4npu2  Jo 
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+Fx$&c osO  Ji(sr) 
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Fx y  sinO  Ji(sr) 
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-r ds 


e-*—\±da 

P 


Displacements  and  stresses 


For  zeroth  order, 
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um(r,z)  —  vm(r,  z)  = 
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For  sin  6  order, 
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A. 2  Dip-slip 


Potentials 


The  potentials  following  the  definition  in  Eq  A.2 
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roo  0.5  sin  25 

a 

3  -0.5  sin  25  4 e"^!2"2*1 

roo  f  —0.5  sin  25— e_Ql2_2*l 


[  -0.5sin25(s2-2/?2)|e-^2-2*l 

00  -0.5  sin  25  (s2  +  02)  £e~al2-2'l 

+  sin  25  s2/?e-0l2-2'l 


Ji(sr)sds 


J2(sr)sds 
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_f71  I  _T7L 

°rz  +  °6z 


rxi 

°rz  ~  °Bz 


4npu)' 


;M0p  f 
Jo 

;M0p  f° 
Jo 


oo  —  sin  26  fsse  “I*  *•! 
+  sin  25 

oo  sin  25  f53c-“l*-*'l 
—  sin  25  fs/?2e-^*-* 


Jz(sr)sds 


Ji(sr)sds 


A. 3  Strike-slip 


Potentials 


The  potentials  following  the  definition  in  Eq  A.l  are 


Mo  iwt  sin  5  sin  25s2  J2(sr) 
47r/)w2  Jo  —  2  cos  5  cos  5f  saJl(sr) 


smusmips.snsrj  -B\z-z  i  s 

e  P|  *  '—as 

—  cos5cos5f/?Jo(sr)  ^ 


Mo  iu)t  y00  sin  5  sin  25s  J^sr) 

- 6  /  € 

47rpw2  Jo  _  CQS  £  CQS  0£pjo(sr} 

Mo  •  y°°  sin  5  cos  25s  Ji(sr) 

- cl  c 

4irpoj2  Jo  +  cos  5  sin  5f/?  J0  (sr) 

-^-re,ut  [  cos  6  cos  6s  Ji(sr)e~^1!~a!,^ds 
iirpu2  Jo  v  '  (J 


Mq  iuit  y°°  sin  5  cos  25s Ji(sr) 
47rpw2  •'O  4-  rns  5  sin  5c/?.7/J.< 


OiiX  C/  WO  AVOUI  I  _  .  |  .s 

V  '  *  ^5 

+  cos  5  sin  5^/3  J0  (sr)  P 


The  potentials  following  the  definition  in  Eq  A.2 


M0  iut  f°°  sin  5  sin  25s2 J2(sr)  i 
47rpu;2  Jo  —  2  cos  5  cos  5fsaJ1(sr) 


47T  pw! 


yoo  cos  i 
^  -si 


oo  cos  5  cos  5s2*3  *2  Ji(sr) 


0  ^  ~  sin5cos25f/?J2(sr) 

roo  T  cos  5sin5f/?^  Ji(sr) 


47rpw 


°_e-<  r  cos  o  sin  tfgp—Ji[sr)  s  ^ 

w2  0  +  sin  5  cos  25  k2J2(sr)  P 
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Displacements  and  stresses 

For  cos  6  order, 

eiut  -oo  —  2cos£o!se~al*~**l 

w  =  - - rM0  /  Jifsrlsds 

4npu2  Jo  +rns£(2*2 


,m  .  ..m 


4npu>*  Jo  +  cos  6 (2s2  -  k2)je-^~^ 

e  »'“*  /*oo  —  2  cos  5  C32e_“llt_z*l 

u  +  v  —  - - rMo  /  J2(sr)sds 

4ttpw2  ;0  +2cos5rs2e-^lz-i*l 


e*wt  -oo  +2cos6fs2e  “lz-2:*l 

=  - - rM0  /  J0(sr)sds 

4W  y0  -2cos5f/32e-^lz-z«l 


°zz  = 


_m  i  _77i 

<7,*  +  <7/9*  = 


c,u,t  z*00  2cos(5^s(2s2  —  A:2)e-alz_z*l 

- - rMo/z  /  Ji(sr)sd5 

47rpw2  Jo  _2cos^f3(2s2-Jfc2)e-^z-z-l 


47rpcj 


r-WoM  f 

•  Jo 


oo  4cos5as2e  “lz  z,l 

—  cos6(y  +  3/?s2)e~^z_z' 


J2(sr)sds 


°°  — 4cos5d£S2e  alz  z*l 


47rpa;2  /o  +  COs<5(2/?3  +  £  +  ^2)e-^lz-z*l 


Jo(sr)sds 


For  sin  20  order, 


wm  = 


4  Trpw2 


oo  sin6fs2e  “lz-z*l 
—  sin  6  fs2e“^z~z*l 


r°o  sin£— e  a'z  z,l 


J2(sr)sds 


um  +  vm  =  - - -M0  /  Q  J3(sr)5ds 

47r^2  Jo  _8in$£e-*l*— I 


«m  -  um 


e*u/t  -oo  —  sin£^e  “I*  *.| 


47rpa; 


;m0/; 


sin5(s2  -2^2)fe-^z~z*l 


azz  — 


_77l  |  ^.771 

<7r*  +  = 


47rpu>: 


47rpw: 


r-Wop  [ 
Jo 

;M0p  [ 

Jo 


<»  —  sin  6(2s2  —  A;2)^e-“lz-z*l 

1  +2sin^/?s2fi-^z-z'l 

oo  —  2  sin  6  £\s3e~“lz_z*l 

+2  sin  6  fs3e-^z~z,l 


Ji(sr)sds 


J2(sr)sds 


Jz(sr)sds 
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Ji(sr)sds 


2sin£fs3e  “I*  *•! 

—2  sin  6  $s(32e~0 1*— *»l 

A. 4  Tensile  Crack 


For  the  purpose  of  short  hand,  the  following  quantities  are  defined  to  be  used 
in  this  section  A. 4. 


Mi  =  M0, 

£i  =  90°, 

II 

5 

62  —  6  +  90°,  and 

Ms  =  ^M0, 

<S» 

II 

O-l 

except  that  Mi  =  —Mo  when  m  =  2,  cos  20  order. 


Potentials 


The  formulation  of  potential  for  tensile  crack  is 

<t>  = 

k=i 

*=1 

3 

'I’e  = 

3 

i’z  =  i 

i=l 

where  the  potentials  with  prime  (/),  which  is  a  single  couple  without  moment, 
are 


Mke 


,iut 


47T/9W5 


r 

Jo 


- v“  - / 

+0.25  cos  2Sk{s2  +  2a2)}  J0(sr) 
—  sin  26k  si11  OcasJ^sr) 
—0.25(1  —  cos  26k )  cos  2  0s2J2(sr) 


.-a  *-*. 


-ds 

a 
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Mkeiut 
4irpu 2 


—0.5  sin  2 6*  sin  0$f3Jo(sr) 

+  {0.25(1  —  cos  26*s) 
—0.25(1  —  cos26k)  cos  26s}  Jo(sr) 


(3 


MkeiuJt  r°°  —0.5  sin  26k  cos  9$(3Jo(sr)  ]  e-p\z-z,\  »  ^ 

Jo  ” 


47rpu;2  Jo  I  -f-0.25(l  —  cos25t)  sin20sJ1(sr) 


Mkeiwt  iwt  0.5(1  +  cos2£jt)$7?*7o(sr) 


47 rpo;2 


-0.5  sin  26k  sin  6sJk(sr) 


e-P\*-‘-\tds 

(3 


Displacements  and  stresses 

The  displacements  and  stresses  can  be  expressed  in  the  form  of  summation  for 


each  orders,  for  short  hand. 


“(r,z)  =  £<(r,») 


For  zeroth  order, 


w?  = 


u?  -v?  = 


{0.25(s2  -  2 a2)f 

gtwt  /•  oo 

47 rpoj2^k  Jo  — 1 0.25  cos  26*(s2  +  2a2)£)  e~“l2~2,l  Jo{sr)sds 

+  (0.25  +  0.75  cos  26k)$s'1e-f,\z-z'\ 

{o.25(«2-2a2)A 

giwt  roo  1  “ 

4nptj2 Mk  Jo  —0.25  cos  26*  (s2  +  2a2) e-Q,k-1!*l  Ji{sr)sds 

+  (0.25  +  0.75  cos  26*)fs/?e-^2_z,l 

f  — 0.25(s2  —  2a2)- 

eiwt  /.oo  *■  “ 

~4npiJ2^k  Jo  -0.25cos2<5*(s2  +  2a2)^}e-al*-**l  J-i{sr)sds 

-(0.25  +  0.75  cos  26*) 
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4ft  pbj' 


_m  I  _m  _ 

°rz,k  "i*  — 


4npu)‘ 


m  m  _ 

°rz}k  ~  Gdz,k  — 


47T  pu‘ 


_£f±^  {0.25(s2  -  2a2) 

f  00 

;Mkp  /  -0.25  cos  26*(s2  +  2a2)f}  e~“lz-z*l  J0(sr)sds 

J  o 

-(0.5  +  1.5  cos  26k)/3s2e~P\z~z,\ 

{— 0.5(s2  -  2a2)fs 

TOO 

;Mkp  +0.5cos2£*(s2  +  2a2)fs}  e~“lz~z,I  Ji(sr)sds 

-(0.25  +  0.75  cos26*)fs(s2  +  0a)e-fl—'»l 

{0.5(s2  —  2a2)fs 

r  oo 

:Mkp  — 0.5  cos  26*(s2  +  2a2)fs}  e~“lz-z*l  J_i(sr)sds 

J  0 

+  (0.25  +  0.75  cos  26k)^s(s2  +  /?2)C-^-z*l 


For  sin  0  order, 


<  = 


«r  +  ^r  = 


4npuj 


;Mk  f 
Jo 


°°  as  sin  26k e  “I*  *•! 


« / ^ 

U*  Vk 


47rpw2 


[  -0.5|(/?2  +  52)e-^ 

/•oo  $\s2  sin26*e-a,lz-z,l 
Mk  / 

— fs2sin26*  e~^z_z* 


Ji(sr)sds 


Ji(sr)sds 


eiut  -oo  —2  sin  26*  cs2e  “lz  z*l 

- - /  J0(sr)sds 

4ttpu;2  Jo  [  +2  sin  24 

e*“*  -.oo  — 2sin26*  £s(s2  + /?2)e_alz~z#l 

- zMkp  I 

4nP w  ;°  +2 sin  2(5*  <rs(s2  +  /?2)e-0l*-*.l 


Ji(sr)sds 


e  r 

<*+<*  =  W^/o 


oo  — 2sin26*  as2e  alz  z*l 


47rpa;2  Jo  +0.5sin 28k  f  (3/32  +  s2)e-0l'-*.l 


J2(sr)sds 


to  to 

°Tz,k  ~  a6z,k 


4npoj‘ 


;Mkp  f 

Jo 


oo  2sin2(5*  as2e““lz-z*l 

-  sin  26k  (0s  +  0.5^  +  0.5/?s2)e~/Jlz-z'l 


Jo(sr)sds 


For  cos  20  order, 


wY1  = 


4irpoj 


;Mk  f 
Jo 


oo  0.25fs2e-alz_z*l 
-0.25fs2e-^lz"z*l 


(1  —  cos  2(5*)  J2(sr)sds 
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«r+«r 


m  _  m 
uk  vk 


4npu: 


4irpu}' 


roo  0.25— 

;Mk  “  (l-cos2  6k)Jz(sr)sds 

’  Jo  _n  okbIc-P  *-*« 


-0.25Ve_/Jl*-**l 

P 


roo  —O.25~€“a‘2r“a'0l 

;Mjt  /  “  (l  —  cos2<5*) Ji(sr)sds 

Jo  -0.25(s2  -  2/32)ie-^-''l  J 


«•“*  /-oo  -0.25(52  +  /?2)4c-“|z_'e'1 

- - -  Mkp,  (1  —  cos26k)J2(sr)sds 

4wpuj*  Jo  +o.552^-^I*-*'I 


_m  ,  m 
ar^,Jk  "T  v$Mtk 


_T7i  _ra 

—  °6z,k 


4‘KpU) 


roo  -0.5fs3e-“li-2'l 

t Mkn  /  (1  -  cos28k)Jz{sr)sds 

Jo  +0.5fsse-^-**l 


4tt  pui 


;Mkp,  I 

Jo 


oo  0.5fsse-alJ'_*'l 
-0.5fs/?V^*-‘*l 


Ji(sr)sds 


A. 5  Explosive  Source 


Potentials 


M0eiut 
4'KpOJ 2 
0 


ah2e~a\z-^-ds 
Jo  a 


Displacements  and  stresses 

Since  the  explosive  source  is  omnidirectional,  there  exists  only  zeroth  Fourier 
order. 


Mne,ut  r°°  „ 

wm  =  — - $h2J0(sr)sds 


um  +  vm  = 


4irpw2 
M0eiut 
4tt  pbJ2 


roo  c 

/  h2—Jisds 
Jo  a 
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_ 

M0eiwt  f°°u2s  T  ,  N  J 
.  ,  /  h  J_i(sr)sds 

input*  Jo  a 

l L 

U  — 

__77l 

°zz  — 

am 

rz 

+  °Bz  = 

Mnpiut  r°o 

- - p  /  2 sh?$JAsr)sds 

4npu*  Jo 

—  am  — 

Mnpiwt  r°° 

4^"/o 

°rz 

°6z  ~ 

It  is  noted  that  the  seismic  moment  Mq  has  been  used  for  source  strength  of 
explosive  source  in  a  solid  medium.  However,  the  quantity,  ^r,  can  be  readily 
converted  into  the  volume  displacement  using  the  relation  h?  =  which  is 
widely  used  in  the  acoustics  community. 
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